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The principal objective of this study is to develop 
a sequence so that it v/ill be possible to analj^se matter as 
available in na,ture, having certain geometric properties 
and environmental stipulations under suitable approximation 
that are not over restrictive-. With this idea in mind. 
Chapter 1 has been devoted for the formulation of the gene- 
ral equations in continuum mechanics. In this chapter, Sec- 
tion 1 deals with the reference frame, characterising the 
variant state of a body, while Section 2 gives the principle 
of material frame indifference, which a physical law should 
always satisfy. In Section 3 the type of forces and 
stresses used in this sequel, has been defined through the 
concept of power generation having conjugates as velocities 
and velocity gradients. In subsequent sections, the govern- 
ing equations for a continuum has been derived through the 
principles of thermodynamics and material frame indifference 
for dipolar stress field wdiich ha.s been summarized in the 
last section. 



XX 


In Chapter 2, the equations "based on finite element 
approximation and parametric differentiation "have been 
developed. The first few sections of this chapter deals 
Y/ith the historical review of finite element technique, re- 
quirements in the formulation and the derived equations 
based on general interpolation formula. These equations, 
in general, will be highly complex and nonlinear with res- 
pect to chosen field variables. To solve this set of equa- 
tions, parametric differentiation technique has been adopted. 
This will reduce the set of equations into another linear 
set of first order differential equations with variable 
coefficients in terms of the incremental field variables 
with respect to a certain chosen parameter. These equations 
together with a standard quadrature technique can be employed 
to solve the problem. 

In Chapter 5} "fche basic equations derived in Chapter 
2, have been simplified for two dimensiona.1 plane problem 
choosing rectangular finite element model. A few problems 
concerning elastic, viscoelastic and thermal cases have been 
solved for numerica,! example. Chapter 4 has been devoted to 
the simplification in the case of axisymmetric solids of revo- 
lution, discretised by ring shaped triangular finite element 
model having only linear distribution of displacements, 
expanded in fourier series in the circumferential directions. 
Examples have been provided by the solutions of a few harmonic 


xxi 


loadings for isotropic and anisotropic linear material v/ith 
linear geometry. Nonlinear cases have been shovm for axi- 
symmetric loading and also for viscoelastic material. Chapter 
5 deals with the pure torsion problem v/ith linear and higher 
order (HOT) stress function distribution over a triangular 
element. Numerical examples have been given for different 
cross sections and nonlinear materials. 

In Chapter 6 and 7, mixed triajigular models have been 
developed for plate (in cartesian coordinates) and shell (in 
curvilinear coordinates) assuming linear distribution for all 
the displacements and stresses. In these two chapters, 
Reissner’s variational approach, or more precisely, simpli- 
fied equations of Naghdi’s formulations (with shear correc- 
tion) have been emplo 3 ?’ed to get the equations for the finite 
element. Moreover, error analysis and convergence criteria 
for the plate model has also been studied by application of 
Taylor series expansion and 'classical order of accuracy' 
method. Finally, concluding remarks and comments on future 
reseaxch have been given in Chapter 8. 



GENERAL INTRODUCTION 


Keeping pace with the progress of science, thfe techno- 
logists want to set up and solve their prohlems more and more 
precisely and accurately. For precise formulation and for the 
Solution of problems mathematics become an indispensable tool. 
But it is also of immense importance to remember that the set 
of equations representing a physical system describes only a 
mathematical model which is used to approximate the physical 
reality to a certain degree. It is usually possible to state 
the physical processes in terms of relations between some 
variables, which generally appear in such equations as phy- 
sical parameters and are to be determined experimentally. 

These experiments become meaningful only when they can be 
performed with sufficient accuracy and iduplicability and 
can be analysed with reasonable precision for the state of 
stress and defoimation in the test specimen under suitable 
approximations which are not over restrictive. The study 
concerning the behaviour of such specimen which is essentially 

a form of matter available in nature, has its origin in the 

1 

concept of thermodynamics . 

p *7 

The. science of thermodynamics may be broadly said to 
deal with (1) energy and its transformations and (2) the concept 
of equilibrium and its evolution with particular 2 ?eference to 



2 


systems involving thermal effects. This is merely a conve- 
nient v/ay of classifying all applications of thermodynamics 
and indicates how broad a field it serves. ’Pure’ thermo- 
d 3 mamics is based on tvro fundamental laws. These laws were 
developed and extensively tested in the latter half of nine- 
teenth century and may be said to rest on a very broad foun- 
dation of experience, A third law has been developed during 
the present century’'; but although of a fiuidamental character 
and based on so'und experimental and theoretical grounds, it 
has nothing like the scope of application of the other two 
laws. It is sufficient for the present purpose to note that, 
v;hereas the first two laws led to the definition of new and 
useful functions, the third law, in essence, merely limits 
the value of one of- these functions (entropy). Since it is 
outside the present scope of stu.dy, henceforth no reference 
will be made about the third law of thermodynamics. 

Ill naturallyr occurring changes or processes are 
’irreversible’ in the sense that they tend to prdSeed in a 
certain direction and ne\rer, vjhen.1eft to thems_elveSj_^ to 
reverse, or go in the opposite direction, like a clock they 
tend to ’run dovm’ and cannot ’rewind' themselves. These are 
matters of common experience. STevertheless, the second law 
of thermodynamics is nothing more or less than a generalized 
statement of such common experience. However, many of the 
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changes in nature are reversible, at least partially. The 
important point is that things can he done only at the 
expense of some other system, which itself runs dora. Work 
must he done to reverse these changes, and the work can he 
obtained only hy allowing a spontaneous change to occur 
somewhere else. To state this another way, it is possible 
to undo the result of a given spontaneous change only hy 
allowing another such change and useftfL work can he secured 
only when the tendency of a spontaneous change exists. 

The tendency to change can he recognized as being 
due to certain ’driving forces’. The change may not appear 
to take place because of some resistance that may slow up 
the change as to make it imperceptahle. If the existence 
of the driving force is recognized, then it is possible to 
predict the direction of the change. 

The end point of any spontaneous process is a state 
of equilibrium in ¥7hich forces are balanced and there is no 
further tendency to 'change. It is important to make a dis- 
tinction between a true, or stable, equilibrium with its 
balanced forces and an apparent, or false, equilibrium in 
which the tendency to change still exists but, owing to a 
high resistance, the rate of change is so small that to 
all intents and purposes no change is taking place, i,e,, 
cessation of motion does not necessarily indicate a true 
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equilibrium. Thus, it is seen that there is difference between 
’true equilibrium’ in thermodynamic point of view and ’equi- 
librium’ from mechanical concept. 

In mechanics, there are two other laws characterizing 
the mechanical state of a system. They are known as Cauchy’s 
lav/s of motion. One of them states the conservation of mo- 
mentum and the other furnishes the conservation of moment of 
momentum, 

Unfortimately, the thermodynamic principles, in them- 
selves, are not sufficient to prove the Cauchy's laws of 
motion. Although from a special theiTnod3rnamic consideration, 
the first lav'^ can be derived, but it is not possible to 
derive for the second law. 'The equations of mechanics des- 
cribe a wider range of phenomena than do the principles of 
thermodynamics. Ho one v/ill contest the balance of mass, 
momentum and energy?', but the existence of a caloric equation 

-j 

of state is an assunption of a more special kind’ .In 
fact the principal objection regarding thermodynamic laws 
lies in the concept and existence of entropy in nonequili- 
brium state. 

As long as a, system is in equilibrium, the classical 
thermodynamics (not statistical) are logical and self con- 
tained. But the situation is essentially different for the 
systems which represent a more wider irreversible nonequili- 
brium process. Here, the question arises whether to give up 
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entirely the concept of entropy and adopt an ’entropy free’ 
til ermo dynamics or to extend the definition of entropy to 
states far from equilihrinm^’ However, even without intrud- 
ing deep into the fundamental validity of thermodynamics, it 
ma3'' he stated that some other hypothesis is necessary to 
derive the Cau.chy’s laws from the principles of thermody- 
namics. The hypothesis adopted here is the objectivity 
principle of physical laws. 

The objectivity principle plays an important role in 
the theory of relativity. Here this principle is applied on 
space-time continuum. In classical mechanics, the velocities 
of different systems are considered to be very very small 
when compared with the velocity of light. In this situation, 
time measured in different intent ial frames is invariant. It 
v/ill seen that when this simplified objectivity principle is 
applied on the laws of thermodynamics, the Cauchy’s laws and 
thermal equation of state will be naturally separated so as 
to ensure the objectivity of the laws. The subjects upto the 
derivation of the general laws of motion and heat equation 
will be dealt in Chapter 1, 

Once these general equations are obtained, it is a 
matter of technique to get a. solution out of them for a 
definite problem. Theoretically it is possible but the 
amount of complexities encountered for the exact solution are 
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virtually prohibitive. Even the classical approximate pro- 
cesses such as as3niiptotic expansions and weighted residual 
techniques can be applied to relatively simpler problems. 

Due to the advent of high speed computers in last few decades, 
the approximate techniques such Ritz, finite difference, 
v/eighted residue, invariant imbedding etc. again have come 
to limelight. But each one of them have got their inherent 
shortcomings. Perhaps, the finite element method may be 
considered as one of the most povi^erfifL and versatile dis- 
cretization presently available among the approximate tech- 
niques. 

Numerical methods, as such, are very useful because, 
when properly applied, they can be used to solve within 
engineering accuracy for the unknovn field variables, of ex- 
tremely complex systems. Historically, difficuJ.t engineering 
problems ha,ve been approached by simplifying the system to 
one that could be solved using either simple formulas or the 
closed form solution of the governing equations for the much 
modified system. Of course, analytical complications may 
arise in many problems. The difficulties, generally encoun- 
tered, are the following: first the overall character of 
the system has been altered, and second, in the zones where 
steep gradients of the field variables exist, the fidelity 
of the solution is questionable. Hence, an accurate model- 
ling of these . situations is very desirable. As the functions 



7 


which the system performs becomes more critical, the need for 
a less approximate method of analysis becomes more imperative. 
Therefore nimerical analysis techniques such as finite ele- 
ment method have developed in an attempt to model the real 
system more closely. 

Using the finite element method, the system is ima- 
gined a,s being divided into a number of finite regions, or 
elements. The behavior of the individual elem-ent is assumed 
approximately. Then the behavior of entire system is assemb- 
led from tha.t assigned to the individual elements. Satisfac- 
tion of complicated geometry and boundary conditions do not 
pose a very great problem. Neither the finite difference 
nor lumped parameter methods possess this advantage to the 
same degree as the finite element method. Even the zone 
of high gradients can be more simply modeled using finite 
elements, though it alv/ays requires special attention. Thus, 
the finite element technique is selected over other discrete 
methods of analysis. 

There are various alternative ?7ays for obtaining tbe 
equation resulting from finite element discretization. By 
far the most popular method is based on variational princi- 
ples W'ith respect ,to minimization of potential or complemen- 
tary energy or Reissner’s energy formulation. However, any 
other variational approach of mathematical physics is equalljj" 
applicable. In this sequel, Galerkin’s technique has been 
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applied in the first few chapters where the approximate equa- 
tions have heen obtained from the general derivations of 
Chapter 1. A part of the torsion anal 3 /'sis vhich has been 
dealt in Chapter 5 , the principle of mmnimum c om.pl ement ary 
energy?' has been utilised whereas the equations for plates 
and shells hove heen derived through Reissner's approach. 

Finite element technique annihilates the space depen- 
dence and expresses the field equations in terms of the field 
variables at discrete points, called nodes. The res-ulting 
equations, in general, will be nonlinear and simultaneous 
which obviate the use of any direct method of solution. Hence, 
necessity arises for the adoptation of a systematic incremen- 
tal approach. A paxametric differentiation scheme has been 
chosen here as the method of stepwise solution of the equa- 
tions. This is adopted because it appears to be quite logical 
and straight forward for application. In a nut-shell, the 
method consists of differentiation of the governing equations 
Vijith respect to a parameter ivhich is assimed to be the only 
independent variable of the solution space where a solution 
is represented by a point in this space. The external pro- 
portional loading in the case of time independent system or 
time in the case of time dependent system may be selected as 
this parameter. The resulting set of equations will he linear 
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and simultaneous with respect to the differentials and can 
he solved. Nov/ it is only left to calculate the actual 
values of the field variables for a small increment of the 
parameter using an 3 / guadrature formula. This way, the solu- 
tion will march out upto a required value of the parameter. 
The convergency of this process is somewhat dubious. But 
the question seems to be related to the order of the quadra- 
ture formula and the step length, both of which can be 
improved in the computer programming at the expense of 
execution time. The finite element and the associated 
parametric differentiation scheme have been presented in 
Chapter 2. 

Different applications to structural problems have 
been shown in Chapters 3 to 7 v/ith numerical examples and 
the general conclusion with comments on future research 
has been presented in Chapter 8. A note for the computer 
programs used for solving different problems has been given 
in Appendix A. 
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CHIPTER 1 


GOVEMING EQUIIIOES 


1.0 Intro (11101:1011; 

Matter is commonly available in nature in the form 
of materials. Inalytical mechanics elude this fact, intro- 
ducing concepts of mass point and rigid body which is 
independent of its constitution. Modern physics concerns 
with the elementary particles of matter. It evades the 
question how a specimen, made up of such particles, will 
behave in natural circumstances. The first approach deals 
with the study of rigid body, whereas, the second approach 
concerns with atomic physics. A theory aiming to bridge 
this gap and to furnish the matter as available in nature 
must describe their mechanical behaviour, deformability 
and represent the definite laws thei' obey under certain 
approximations which are not over restrictive. 

The complete mathematical statement describing the 
state of a deformable medium consists of momentum and mass 
conservation, thermodynamic lav/s, certain geometric pro- 
perties and environmental stipulations expressed as boundary 
and initial conditions on loads, displacements and/or heat 
flow. The variables commonly used to describe these changes 
of state are stress and displacement components, temperature 



12 


entropy and heat flux. The variant state of a "body must he 
specified through frames of reference, which has been dis- 
cussed in Section 1. Moreover, the laws of nature should 
not he special for a frame (or an observer). This objecti- 
vity principle is dealt in Section 2, In Section 5, different 
types of forces relevant for subsequent analysis have been 
defined through the principle of power generation. 

The whole domain of classical mechanics, which is 
the basic background of the present v/ork, is based on the 
approximate ne’rtonian concept and the classical law of non- 
interconvertibility of mass and energy. But there exists 
various forms of energy. They are amenable to change their 
form during a change of state or so to say, during a process, 
iny process will be governed by thermodynamic laws (by defi- 
nition). In actuality, as any other laws, they will put 
restrictions on the inter diffusion of various form of 
energy in matter, e.g. , mechanical, heat, electrical mag- 
netic etc. In this sequel, only the first two types will 
be considered. 

The ultimate idea of continuum mechanics is to pre- 
dict the behaviour of ma,tter during a process, as closely 
as human knowledge permits. Since the present progress 
speaks in all certitude that matter is discrete and energy 
transformations for the elementary particles follov/ only 
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probabilistic certainty even for very slow and gross motion, 
it leads to a direct contradiction of tlie principle of cla- 
ssical thermod 3 maiiiics and consequently of continuum mechanics. 
A happy reconciliation can be made if statistical thermody- 
namics is introduced. But the added complexity due to the . 
introduction of statistical mechanics may not be commensurate 
with the accuracy achieved in a sense of solution of general 
practical problems. However, though the present motivation 
is in no 'na.j so ambitious, a procedure will be adopted so 
that more and more exactness can be incorporated to meet 
the ever increasing accuracy required by future technologists. 

With this idea in mind. Section 4 is devoted for 
thermodynamic concepts and conservation lav^s. In Sections 
5, 6 and 7, the balance of energy will be considered in 
details and the laws of motion ?^rill be derived from the 
principle of material frame indifference. Section 8 deals 
with the equation of thermal state, which will be deduced 
from the expression of entropy. The general physical laws, 
in themselves, do not suffice to determine the deformation 
of a body. Before a determinate problem can be formulated, 
it is usually necessary to specify the material of which the 
Dody is made. In continuum mechanics, such specification is 
stated by constitutive equations. A very brief discussion 
on constitutive equations will be made in Section 9. S'oi’ 
convenience and future references, the important equations 
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v/hich will le needed in subsequent auelysis^ will be 
summarised in Section 10. 

1.1 Kinematic Preliminaries: 

To study the motion of a continuum, it is necessary 
to establish the correlations between different phases of 
the motion which the continuum undergoes with respect to 
time and relative to a reference state. The description 
of such a reference' state must be characterized through a 
set of time dependent coordinate frames. 

let B be a three dimensional manifold of material 
points within the material volume v and bounded by the sur- 
face s at initial time t = t . The rosition of a material 

o 

point X in this region may be denoted by one of the coordi- 
nate frames, i.e., a rectangLilar cartesian coordinate a 
curvilinear coordinate X^ or by & position vector X, where 

X^ = X^ (z\ Z^, Z^) = X^(Z^) 

X = X (X^) = X (X- (Z^)) 

and the indeces having latin miniscules (i, j, k etc.) take 
the value of 1 , 2 and 3 and refer to the reference coordinate 
system, in the undeformed state. 

A motion of the body B is the family of mappings with 
respect to time parameter t of B into Euclidian space B.^., 
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The image, 

X = X (X^, t) (1.1-1) 

of X is called the position vector of X at time t. In com- 
ponent form, Eq. can he written as, 

x^ = x^ (X^, t) = x^ (Z^, t) (1.1-2) 

having four independent variables t and X^ or Z^, i = 1,2 
and 3. 


It is assumed that mappings (1.1-2) are single valued 
and have continuous partial derivatives with respect to their 
arguments to whatever order desired except possibly at isolated 
singularities like points, curves and sixrfaces which necessi- 
tate special attentions. Eurther each member of (1.1-2) has 
unique inverse of the other in a neighbourhood of a material 
point X. Hence position vectors X and x may also be expressed 
as function of x^ at any instant of 'time. But to save space 
and to avoid any imbroglio ■v'.tiich may appear in subsequent 
derivations for finite elem^ents, all quantities, in this 
sequel, will referred to undeformed state X^ as has been 
shoTO in Ea^s. (1.1-1) and (1,1-2). Further axioms of indes- 
tructibility and impenetrabilitj'- of matter will lead to the 
inequality, 



0 


(1.1-3) 
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iny deformation of continutun, treated in this variant, will 
te tacitly asstraied to respect this inequality. 

The covariant components of base vectors and metric 
tensor for the undeformed and deformed configurations may 


be d 

efined 

as, 






Gi 

= \ (x3) = 

SX 

8X^ 


( 1 . 

1-4) 


Si 

= ii (X^, t) 


8x 

ax^ 

( 1 . 

1-5) 

and 

G . . 

ID 

= G. . G. 



( 1 . 

1-6) 


Siq 

= gi . g^^ 



( 1 . 

1-7) 

The 

contravariant components 

of the corresponding 

quantities 

can 

be obt£ 

lined by solving th 

e equations 





• = ^"k’ 

G^ = 

^ G^^ G. 

0 

(1. 

1-8) 

and 

p-iO 

o 

* S3X " 

g^ = 

ii 

0 ej 

(1. 

1-9) 


where 5^^ denotes Krone cker delta. 

Instea.d of x, the position vector of time t can also 
be expressed in terms of the deformation vector U, such that 

X = X + U (1.1-10) 

where TJ = U (x^, t) 

and tJ (X-, 0) = 0 (1.1-11), 

i.e., time is reckoned from the undeformed state. The 
covariant and contravariant components of the deformation 
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vector are, 

U = (1.1-12) 

and the classical strain tensor can he defined as, 


E. . = 

10 


_L 

2 

1 


^®i3 ^'ij^ 


k 


( 1 . 1 - 13 ) 




.) 


V-liere, har denotes cova,riant differentia.tion with respect 


to metric G. .. 

10 

It is worth mentioning here, tha.t strain plays the role 

of an intermediary to correlate through thermodynamics, the 

auctions and their corresponding response within matter. It 

will he seen later that certain quantities are necessary as 

arguments of functions so as to satisfy certain basic postu- 

1 20 

lates of physics. Ihere are a host of quantities ’ , vrhich 

might he acceptahle and, as such, they may he called strains, 
hut their fundamental dependence on displacement gradients 
are not important except to satisfy the physical aspects for 
which it is meant and relation (1,1-13) is only one of them, hut 
?jill he used here to conform to traditions. 
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1.2 Ma^terial Frame Indifference: 


It has heen seen that to specify’- any event in a system, 
it is necessary that the event has to he located with respect 
to a reference frame providing information for space as well 
as time. The specification of a reference frame is not miique 
and depends upon different observers. But the fundamental 
mea.surable quantities, such as distances, angles and time in- 
tervals are independent of the observer. Consequently, any 
change of frame must preserve them along with the temporal 
orders of all the events. The most general relation for such 
a change of frame may be expressed as, 


X* (X^, t*) = R(t) X (X^, t) + D(t) 
t* = t - a 


( 1 . 2 - 1 ) 


where (x , t ) represents an event in one frame and (x, t) 
in another. D(t) is a constant vector representing a point, 
R(t) is a proper orthogonal tensor and a " iS' a real number. 
D(t) and R.(t) ma.j?' be assumed as smooth fimction of t. It is 
not difficult to visiialize that the two frames are related 
by a rigid transformation combined with a time shift and 
differ only by an Euclidian displacement. In caxtesian co- 
ordinate system, Eq. (1.2-1) may be expressed as. 


*1 


* 


(z-% t"-) 


R\(t) x^(Z^^, t) + D^(t) 


( 1 . 2 - 2 ) 


t = t - a 
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Two motions of a given medium related by tbe Eq. (1.2-1) or 
(1.2-2) are said to be equivalent. 


Eor a change of frame, transformation laws for scalar, 
vector and tensor are given as follows: 


(a) 

(•b) 


or 

(c) 


or 


A scalar remains unchanged. 

A vector is transformed into, 

V* (X^, t'') = R(t) , v(X^, t) 

Y*^(Z^, t*) = R^.(t) v3(Z^, t) 

0 

A second order tensor is transformed into, 
S*(X^, t*) = R(t).S(X^, t).R^(t) 

S*i.(Z^, t*) = R^^(t) S^^(Z^, t) R.^(t) 


(1.2-5) 


(1.2-4) 


Junctions and fields whose values are scalars, vectors and 
tensors will be called objective or frame indifferent, if 
both dependent and independent variables transfoim according 
to the above laws. 

Equations governing the motion of a continuum will 
be derived primarily from a postulate that a physical law 
will always satisfy the principle of material frame indiffe- 
rence, applied on the balance of energy principle. A short 
bistort’' together with a general discussion on the principle 
of material frame indifference has been given in references 1, 
2 and 3. 
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1,3 Definition of forces: 

A material body -under tbe influence of external and 
internal forces undergoes deformations. On a quantitative 
basis, there may be three types of forces which can influ- 
ence a medium. 

A. Extrinsic forces: 

These forces are due to external effects and regarded 
as acting upon the particles comprising the body. The body 
force due to uniform gravity is a relevant example. 

B. Mutual forces: 

These forces arise within the body and are regarded 
as acting upon pairs of particles having equal magnitude but 
opposite directions. Therefore, the resultant internal force 
is zero, force between two charged particles is an example, 

0. Contact forces: 

These forces result due to action of one body on 
another through the boundary surface, lines or points. Stress 
field and surface loads belong to this category. 

The quantitative definition of forces can be derived 
easily from the conce-pt of worh rate. However, it is to be 
noted that the definition of forces is not unique nor straight 
forvmrd* The basic idea is that certain affinities (e,g. 
velocities being acted upon by corresponding fluxes (forces, 
P ) produce po?/er (i.e., rate of work W = P fj_). Hov/, 
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conversely if is a contravariant vector and is an arbi- 
trary covariant components of a set of velocity and if the 
scalar P^V^ is a rate of work corresponding to the velocity 
vector, then P^ may he defined as a force in the direction 
of Thus defined, force may be assmned as a function of 

position vector, velocity and time. In continuum mechanics, 
this idea may he easily extended to define forces of more 
general nature 


Surface Porce: 

let an independent set of kinematic variables with 
respect to motion (1.1-2) he denoted by a tensor function 

X. . . under coordinate transformation, where, 

111 ••• 


^ii 5 = ^ii. i ’ -fc) » P = 1,2,3.. (1.3-1) 

In addition, if tensors (1.3-1) respect the principle of 
material frame indifference, then they may he designed as 
simple 2^ pole displacement field or in general multipolar 
displacement field"^. How, consider a surface s whose unit 


normal at point X is n^. If P 




IS a 


tensor of order a + p + 1 and if for all arbitrary 2^ pole 


velocity gra.dients x. . 


1* • --"pi '^1'* 

p • 


, the scalar 




X 


i.i^* • *^8' ^1* * *^oc 


is a rate of work per unit area of s, then tensor 
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P 


• # f lo* ^ 


-j • • * 3 „ 0'.+ B 

may be called a surface force 2 ' -pole 


of the (p + 1) th kind per unit area of that surface^. When 


0, P 


i;0 


1 * * 


a 


denotes the surfs,ce force 2*^-pole of the 


first kind and when a = 0, P ^ is called a simple sur- 

face force 2^-pole of (p + 1) th* kind, a = p = 0 corresponds 
to the classical surface force P^. Various other degenerated 
cases can be obtained depending on the order of the tensor 
X. . . 5 order of gradients a and the special properties 

Xi1 ... ip 


imparted to x. . 


11 ^ . . . 1 


P 


(for exiample = -x^. Such 


1 


rii- 


special cases and their relations with the director theories 
has been discussed by Green and Naghdi^. In this study, only 
the dipolar (a = 1) force of first kind (p = O) will be con- 
sidered. Hence, the total rate of work over the surface s 
by the monopolar surface force P^ and dipolar surface force 
p^^ corresponding to the velocity V^. and velocity gradient 

3 ji 


'^’-=1 . will be given by, 




( 1 . 3 - 2 ) 


Stresses: 

Again consider a surface normal to X^-aocis in 
the Tnitial configuration. The corresponding contravariant 
components of the stress (monopolar) and hyperstress (dipolar) 
of first kind may be respectively denoted by T^^ and 
if the Corresponding rate of work done is given by, 

T^^ V. and T^^k Y (1.3-t3) 

J 1 J - 
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per unit area of s. for any arbitrary velocity V. and velo- 
city gradient "^j first index i of the stress and 
hyperstress components indicates the surface on which 
they act. 

Body forces: 

If and are contravariant tensors and and 

Y^ j i arbitrary set of velocity and velocity gradient and 

if the scalars P^V. and P^^V.i . are rates of work per unit 
- 1 _ 3i 1 

mass, then P^ and P^^ are respectively called body forces 
of monopole and dipole of first kind per unit mass. The 
total rate of work of the body forces distributed throughout 
the volume v is given by, 

f (pH. + dv (1.3-4) 

where is the density of the mass at time t = t^. 

1,4 Thermodynamic Preliminaries: 

The basic idea of continuum mechanics is to formu- 
late theories which will reflect, during any change of state, 
the overall expositions of the microscopic behaviour of 
matter. Moreover, for obtaining the complete system of eq_ua- 
tions describing the motion of the medium, it is necessary 
to know the general laws governing the changing phases of any 
medium and the characteristics of the particular medium. The 
necessary informations should come from the general principles 
of thermodynamics. 
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A thermo dynamic system may he defined as a system 
in nature whose action on its surroundings proceeds hy way 
of output or absorption of work and heat. This definition 
applies even to a nonequilibrium process for which the equi- 
librium state is only a particular case. It characterizes 
thermod3niamics by operations, i.e., the exchange of work 
and heat between the systems and its surroundings and it 
relates these quantities with the reaction of the system, 
which are deformations, the temperature and their gradients. 
Now a system may include discretely the fundamental particles 
of matter or it may specify larger aggregates exposing only 
the average reactions of particles which is the basis for a 
mathem.atical model of a continuum. 

Since this study is based on continuum mechanics 
approach having no chemical and electromagnetic effect on 
the medium, it will be assumed that any physically small 
fundamental elem.ent of the medium can be reganded as a ther- 
modynamic system, kor each of these small elements, the 
mechanical concept of internal states are defined and each 
successive state of evolution is characterized by a finite 
number of defining parameters. The exchange of work and 
heat will be not only over the surface of total system but 
also with the surroundings of any specified subsystem. 

Any fundamental element at any stage of evolution 

5 

will follow three basic laws'^; 
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1 . Conservation of mass: 

The mathematical expression for conservation of 
mass may he expressed as 

. k 

P (x, t) = Po(X, t )/i^ I (1.4-1) 

where p (x, t) is the density of mass at time t for a point 

X which was at X at time t . It will he assimed here that 

0 

any deformation from the reference state will always satisfy 
Eq. (1.4-1) 5 or in other way, this relation may he accepted 
as a definition of density at time t for the point x. 

2. Balance of Energy: 

Balance of energy may he expressed as, 

K + i) = P + ZE^ (1.4-2) 

a 

where X, E and P represent respectively the kinetic energy, 
internal energy and the power supplied hy the external mecha- 
nical forces. represents the mechanical equivalent of the 

ath kind of nonmechanical energy per unit time, e.g. , heat 
energy. The dot denotes differentiation with respect to 
time. This law may he designated as the first law of ther- 
modynamics. 

In a given physical prohlem, P, and K are gene- 
rally well defined. The remaining quantity, i.e., the inter- 
nal energy may he looked upon as the quantity which halanees 
this equation. Since internal energy is also an additive 
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function, in terms of the internal energy density e per imit 
mass, it can he written as, 

E = / p e dv (1.4-3) 

V 

The internal energy density e is a material frame indifferent 
scalar and depends only on thermodynamic variable. 


3. Second law of Thermodynamics: 

If M he the number of mechanical parameters 
which influence the state of motion of a material element 
having specific internal energy e, then the second law of 
thermodynamics asserts that for any given values of e and 
^ ^ in a certain domain, there exists a quantity called 
entropy whose specific value n depends on e and v i.e,, 

n = n (e, V x^) (1.4-4) 

The differential form of Eq. (1.4_4), 

dn = x ^ a (1.4-5) 

a a 

gives the concept of 9, called the absolute temperature, 
though the equation, 


1 

e 

The two functions of 
perties : 


§1 (1.4-6) 

state 0 and n have the following pro- 


I. 9 is a positive number which is a function of the 

empirical thermometric temperatu 2 ?e. 
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II. The entropy is an additive function and the sum of 
entropies of the parts is equal to the entropy of the whole 
system. 

III, Change of entropy can he affected by change of heat 
energy and by production of entropy inside the system. If 

•¥r 

dQ denotes the increase of heat energy at temperature 0, 
then the increase of entropy can be written as, 

* 

dn = + (1.4-7) 

where dn^^ is the production of entropy inside the system, 

17. Por any process in nature the production of entropy 
inside a system can not be negative, i.e., 

dn ^^>0 ( 1 . 4 - 8 ) 

The equality v/ill correspond to a reversible system while 
the inequality stands for irreversible process. 

7. The entropy function must fulfil the Glandsdroff- 

O 

Prigogine theorem , 



'Sow if it is assumed that all thermodynamic functional 
relations are differentiable as many times as needed and are 
invertible to yield one variable as a function of the others, 
then from (1,4_4), e can be written as, 
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TTae physical dimension of v ^ axe that of mechanical xmits 
(say displacement gradients and higher gradients), but other- 
wise arbitrary. The temperature 0 and the thermodynamic ten- 
sion T®- may be defined from (1.4-10) by, 

e = 5e ^ p_ (1.4-11) 

a 

It will be assumed that this temperature is the conventional 
temperature. The thermodynamic tensions will be shown later 
to be related with conventional stresses. The differential 
form of the Eq. (1.4-10) can be written as, 

¥ 

de = edn + 2 x® a.v (1.4-12) 

j a 

a=1 

Since this equation is valid for any arbitrary differential 
increment of n and v^, as a special case, it will be valid 
for the actual motion of the particle X. Hence, the corres- 
ponding rate equation can be written as 

M 

e = en + 2 x“ V (1.4-15) 

' — -I 06 

a=1 

Sometimes, instead of internal energy density e, it is 
more convenient to use Helmholtz free energy per unit mass, 

^ , which may be defined as, 



( 1 . 4 - 14 ) 


( 1 . 4 . 15 ) 



29 


v/here, 

= ’I' (e, V x^) 

.The name free energy is used because, it is that portion of 
the energy which is available for doing v;ork at constant 
temperature, 

1.5 Energy Balance: 

Consider an arbitrary material volume v of the con- 
tinuum boimded by a surface s having unit outward normal n 
in the reference configuration. From Eq^. (1.3-2) and (1.3-4), 
the rate of work done by the surface forces over the surface 
s and by the body forces within the volume v is given by, 

P = / (P^T. + P^3 Y ) as 

s ^ > 

+ / (P^Y^ + Y^I^) dv (1.5-1a) 

how if the dipolar kinetic forces, i.e., inertia forces due 
to gradient of velocities are neglected, accepting the kine- 
tic pov/er K in the conventional form, then, 

i = / p (Y^ V ) dv (1.5-1h) 

V ° 

Restricting the derivation for the situation where 
nonmechanical power supply will be only from thermal consi- 
derations and there will be no multipolar heat flux, the 
supply of thermal energy per unit time by heat flux h across 
the surface and by distributed energy sources in the body 
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per unit mass q, is given by, 

2 P = / p q dv + / Ms (1.5-1c) 

a V s 

Heat energy b will be absorbed by tbe material and will be 

supplied by radiation from tbe external system across the 

surface s. Farther, let 6 be a local absolute temperature 

and let be tbe conventional beat flux vector through the 

surface s into the volume v. Substituting expressions for 

K, P and E F in tbe balance equation (1.4-2) and after some 
a 

rearrangements of terms, it yields, 

/ P- (vH + e) dv = 

+/ (h + pV. + P^^V.j .) ds 

In the next step, relations between the surface forces and 

stresses together with tbe beat flux h and heat flux vector 

over the surface s is established. The method adopted 

4 5 

here is similar to that given in references 

Consider the volume v in Eq. (1.5-1) to be an internal 
tetrahedral element bounded by a plane with arbitrary unit 
normal n^ and by planes through the point X parallel to the 
coordinate planes. This element would be subjected to the 
internal stresses and well as heat flux vector 

'and external surface loads P^,' P^3'“ §pg-~^;g=fe5n^'~beat flux 
h. If ds is the area of tbe plane of the tetrahedron normal 


/ P„(q + 

( 1 . 5 - 1 ) 



to and ds^ is the element of area of the plane of the 
tetrahedron normal to then, 


dSf = n^ds ‘ (1,5-2) 

If the volume of the tetrahedron shrinks to zero preserving 
the orienta^tions of its planes, then the volume integrals in 
(1.5-1) vanish and hence this equation can now he modi- 
fied to account for the internal stresses and heat flux vec- 
tor to the form, 

(P^ - n. T^^) V, + (p3^ - V, > . 

J X K 1 j 3 

+ h - n^ = 0 (1.5-5) 

For an arbitrary constant rigid body velocity, stress tensors 
1^^, surface loads P^, P^^, heat flux h and heat flux 

vector remain unaltered. Hence, 

(P^ - n. Y = 0 (1.5-4) 

where is now the arbitrary rigid body velocity. Hence, 

pi = n. 

D 


and Eq, (1.5-3) becomes. 


(p3i _ n^ T^^^) + h - n^Q^ = 0 (1^^-5) 


Again Y ^ j ^ can be written as 

Y. I ■ = A. . + . 

i|3 X2 10 

where the symmetric tensor (deformation rate) 


(1.5-6) 
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and antisymmetric tensor (spinor), 

=' -^3|iV2 (1.5-8) 

Substituting (1.5-6) in the Equation (1.5-5), 

(P^i - nj^ + h-n.Q^ = 0 (1.5-9) 

hext supplying a rigid body angular velocity and observing that 
qI ^ g_j,Q unaltered but W. . changes to 

X j X 3 

W. . W. . + 2^ , . (1.5-10) 

ID ID ID \ ^ / 

and Eq. (i.5_9) yields, 

(p3^ - nj^ T^^^) = 0 (1.5-11) 

Since Q . is an arbitrary antisymmetric tensor, Eq. (1.5_11) 

J 

can be written as, 

p3^ _ (pkji _ ^ Q (1.5-12) 

and (1.5-9) reduces to, 

(p3^ - n ^ . + ^ _ n = 0 ■ (1.5-13) 

il. J -L 

If P^^ and 1 ^ 3 1 qg split into S 3 ninnetric and antisymmetric 
parts, Eq. (1.5-12) and (1.5-13) can be written in a slightly 
different form. Thus, 

pi Di I _ pkj Di 1 = 0 (1.5-14) 

and (p^^^) _ + h - n^Q^ = 0 (1.5-15) 

where the indices enclosing the braces represent the symme- 
tric pant and within the square brackets, the antisymmetric 
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parts. In deriving Eq. (1.5-15), the S3ni!metry criteria of 
A. . has been taken into consideration. 

X J 

Now consider materials for which 
on the arbitrary boundary surface of the body does not depend 
explicitly on the velocity gradients and hence does not depend 
on A. .. Then, since the element of the latter quantity can be 

X J 

chosen arbitrarily and independently of each other, subjected 
to symmetry restrictions, Eqs. (1.5-14) and (1.5-15) yield, 

(p3^ - n^ =0 (1.5-16) 

i.e., p3i = nj^T^3i 

and h = n^Q^ (1.5-17) 

With the help of Eqs. (i.5_4), (1.5_14) and (1.5-17), 
the balance equation (1.5-1) can be v/ritten as, 

/ p W dv + / p e dv = / p (q + E^Y. + .)dv 

+ f n^Q^ ds + / (n^^T^^ Y^ + n^T^^^ Y.|^) ds 

^ (1.5-18) 

By transforming surface integrals to volume integrals by 
usual way after making appropriate smoothness asstanptions, 

Eq. (1.5-18) can be modified to, 



(1.5-19) 
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After some rearrangement of terms, Eq. (1.5-19) yields, 

\|3 - ->0 "" ^0 

+ 1^33- 1 dv = 0 (1.5-20) 


This balance equation should be invariant under a change of 
rigid body constant velocity. Since due to this change, in 
addition to T^^, pi j^i g^so remain 

invariant by uniform rigid body velocity, Eq, (1.5-20) 
reduces to, 

/ +Po ^i ° (1.5-21) 


and 


/ p„; 4 t = / V.| . + 5l| .) dv 


V 


where, t*^^ = ^ ^o ^ 


( 1 . 5 - 22 ) 

(1.5-23) 


Since this is true for any arbitrary volume, dropping the 
integral, the equation may be written as, 

+ Po^’^ - Po = 0 (1.5-24) 

and p„4 = t*3i Y.|j + V.|j^.p^ 4+ (1.5-25) 

• n 2. 

Again, if additional assumptions are made that e, 
and q remain unaltered by superposed uniform rigid body angular 
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velocity, the body occupying the same position at time t, 
then, 

^ = 0 (1.5-26) 

This shows that the antisymmetric part of t*^^ does not pro- 
duce any work. Thus Eq. (1.5-25) will yield, 

A.J + qjjj. + Pp 1 + (1.5-27) 


Eq. (1,5-24) furnishes the classical equation of motion derived 
from Cauchy's first law , i.e., momentum balance principle and 
Eq. (1.5-26), which being true for any arbitrary reduces 

to , 

(1.5-28) 




provides moment of momentum balance equation, i.e,,. Cauchy’s 
second law of motion. At this stage, no further progress is 
possible unless the independent variables in e, i.e,, are 

assumed explicitly, 

1.6 Internal Energy: 

In classical theory, the elastic energy, which is the 
internal energy for an elastic medium in isothermal condition, 
is given by the form, 

E = / e' dv 

V 

where e* is the energy of deformation per unit of reference 
volume. ITow e’ at any point X and at time t is determined by 
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instantaneous configuration x^(X) of an arbitrary small 
neighbourhood of point N(X) containing X. Suppose that x^(X) 
has p number of derivatives in that neighbourhood. Then the 
relative position vector of the point X and any other point 
X in lj(X) Can be written using Taylor’s series^, in the form, 

x^(x') - z^(X) = (X) dX^ + ^ (X)dx3dX^ 

+ ... +;^x^|. . . (X)dX^''dX^^...dX^^’+0(6^^^) 

P- !^1^2*** p 

where 0(6^"'’'^) represent a term of order p+1 in the diameter 6 
of i'(X), In classical theory the energy term is given by. 


e’ = e’ (x^j^, X) 


( 1 . 6 - 1 ) 


The obvious and natural generalization of (1,6-1) is, 

• /_i -_i _i 


E = / e (x^ 

V 


X ..., X 


. X)dv 


in 

Eollowing Noll’s terminology, if N is the order of highest 
gradient present in the energy density function, the material 
may be termed as grade N material and when N is greater than 
one, it is nonsimple material. In the light of this definition, 
the materials referred in the classical theory and represented 
by the Eq. (1.6-1) are simple materials of grade one. 


The Eq. (1.5-25) which gives the rate of increment of 
internal energy has second derivatives as the highest velocity 
gradient which has actually followed from the assumption that 
only dipolar forces are present in the continuum. It 
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automatically leads to the conclusion that atmost the material 
may be of grade two. Thus in a reverse way, the assumption of 
grade two material will necessitate the presence of dipolar 
forces and any higher order forces cannot produce work and 
hence can be completely disregarded. In this situation, apart 
from the variable n which is a scalar, v should span the 
whole space upto second derivatives of the displacement ftuac- 
tions. Hence, the specified internal energy can be written as, 

e = e (p , x^l^, X) (1.6-2) 

But this expression will have to be restricted so that the 
scalar function e will be invariant under rigid body rota- 
tion. There is no unique way of representing this function 
in invariant form. However, Toupin^ has shown that all these 
conditions will be satisfied if it is assumed that, 


e - e (n , X) 


(1.6-3) 


where E. . is classical strain tensor and 

X J 


E. = E. .1. + E., I . — E 1 . 
lOk 13 jk ik] 3 3kji 


(1.6-4) 


The dependence of the density function e on X in the Eq.( 1.6-3) 
is only to account the nonhomogeneity of the material. How, 
if thermodynamic tensions are designated as, 


and 


- 


6e 


0 


eT 


(k3)i ^ p 


0 




8e 

i j 3k 


(1.6-5) 
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then, the time rate of internal energy can he written as, 

Pc® = Pc ep * hUk 

Since Eq. (1.6-6) has to he invariant under superposed uni- 
form rigid body angular velocity, 

gil i-3 i = 0 (l,6-6a) 

e 

and 




To get a more clear picture regarding the nature of 
the function , eliminating e from the Eqs, (1.5-27) and 
(l.6-6h), the expression for the rate of entropy production 
can he written as. 


p ^ e n = (t 


*(iD) A 

" e ^ ji 


+ (j(kj)l _ ^5,(kj)l) (1.6-7) 

+ Po 4 + 

Since thermodynamic tensions and recoverable , the production 
of entropy, is due to excess of mechanical work over the 
recoverable parts and nonmechanical supply of energy. Desig- 
nating, 


and 


m(i3) — +*(^3) m(f3} 

d^ " ^ ■ e-^ 

$(k3)i == Qi(l5:3)i _ 5(kD)i 


d"" * ^ e' 

Eq. (1.6-7) becomes, 

"^(13) _ + iji(k 3 )i Y 


(1.6-8) 


Po9n = ^T 


3i d 


• I *1 + P 1 + Q^i -t 

Lj^k o ^ ^ ji 

(1.6-9) 
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The prefix 'd* is used to denote that the q^uantities are of 
dissipative nature, i.e., are not recoverable. 

Since entropy n depends upon thermodynamic path, it. 
is possible to get its constitutive equation which will be 
a fimction of temperature 0, mono and dipolan displacement 
gradients and velocit 3 r gradient upto a certain order (with 
respect to time, e.g. , rth order dipolar velocity gradient 
is 6 past history of the thermodynamic 

path. The main point is that the dissipative stresses and 
hyperstresses can be obtained from the constitutive equation 
for n provided the entire past thermodynamic path is known. 

Substituting Eq. (1.6-9) into (l.6-6b), internal 
energy expression can be written as, 

i(k3)i 


where T 


; =>( 13 ) 

do) ^ ^(iD) 


e = T"^^^ ^ + q1 (1.6-1G; 


/ 


d^ 


( 13 ) 


( 1 . 6 - 11 ) 


T(k3)i = 5(k3)i ^ '^(k3)i 

e d 


Comparing Eq. ( 1.6-10) and (1.5-25) and recalling that the 
first and second order velocity gradients are independent 
kinematic variables, it can be concluded that, 


r-h 


*(iD) m(i3) 


T^-^^d A.. = 0 

J 


and 




( 1 . 6 - 12 ) 


and consequently, 

t*Cl3) =, 1(13) ani T(k 3 )l 


i(kj)i 


(1.6-13) 
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to 3?e4uc« the mimher of rariablea^ asmwtaciC 

of the dipolar stresses, juoy -{je replaced, hy 

5 (kj)i (i. 6-15)2> when using Eqs, (1.5-28) and (1.6-15)^ » 

the symmetric and antisymmetric parts of the monopolar stresses 
can be written as, 

mCoi) = J(3i) 5(k(3)i) p(Di) rplk(j|i), 

X ^ }k ■ Po ^ “ [k 

(1.6-14) 

and t 1 3i I = _ jCk |3 )i ^ p^pl ji I _ kj j I i 

(1.6-15) 

In the next section, it will be shown that all the 
unknown variables involved in the governing differential equa- 
tions (1.5-24), (1.6-14) and (1.6-15) are not determinable 
independently, i.e., the set of differential equations do 
not form a determinate set. 

1,7 Indeterminancy of Stresses: 

Apart from temperature e, the total number of unknowns 
in the system is twentyone. This consists of nine nine 

t 1 ^3 li and three displacements U^. The symmetric parts of 
i.e., are no. longer independent unknowns, since 

they have been replaced by from constitutive equations, 

which are functions of displacements. On the other hand, the 
number of governing equations are twelve which are three 
from (1.5-24) and nine from (1.5-14) and (1.5-15). Hence 
total system in mechanical parts is not detenninate. 
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If the expressions for 1 from (1.6-14) 

and (1.6-15) are substituted in the equation of motion 
(1.5-24), the terms 1 j and fl 3 j ^ Ij vanish 

because of the antisymmetry of the indices j and k where as 
the covariant differentiation with respect to j and k is 
independent of order and hence are symmetric. Again work 
rate due to this part of stress, 

is also zero. Although, this part does not contribute to 
the equation of motion, nor produces any work, it plays 
important role for the boundary conditions Por the 

sake of simplicity, this part will be assumed to be totally 
nonexistent within the body and also on the boundary. 


With this -limitation, the governing equations may 
be collected together and written in the integral form as 
follows : 


f 


+ p 

- 

p •) Y. dv = 0 

(1.7-1) 

Y 


0 



f 

Y 

( D i ) 

_ f (ji) 


m(k(o)i) ^ p j,(3i) 

j k 0 

) A . . dv = 0 

^ tJ 






(1.7-2) 

f 

Y 

(fl ' 

I 4 

3) 

^ L + p pi ) W. . 

|k "^0 ^10 

dv = 0 


(1.7-3) 


The corresponding conditions to be satisfied on the boundarie 
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are , 

/ (P^ - n- Y. ds = 0 

s 3 ' 1 

/ (P^3i) _ n, = q (1.7-4) 

g 1C 13 

and / (p! ! - n, I 3 ) i 1^ 77_ ,3^g = q 

s Jc 13 

Combining Eq, (i.7_i) with (1.7-4)^, the overall 
momentum balance over an arbitrary small domain can be 
written as, 

f _ p^yi) Y^ dv 

+ / (p^ _ n. Y ds = 0 (1.7-5) 

s 3 1 

How, since Eq. (1.7-2) and (1.7-3) do not contribute any 

new develotaoent , it is aidirantageous to cojablne them toother 

with (1.7-4) 2 and (1.7_4)^ to yield, 

f (iJi . i(n) J(k3)l + p pjl) V., . dT 
V jk 0 ^ il3 

+ / (p3i _ _ ds = 0 (1.7-6) 

In these equations, E^ and are body forces, are 

monopolar unknown stresses, P^ and P^^ are surface forces 
and are unknown velocities, ill the above quantities 
are functions of space and time. Moreover, 1^3^^ and T^^3)i 
are thermodynamic stresses supplied from constitutive equa- 
tions. 
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1.8 Equation for TEeCTal State: 

The equation for thermal state can be easily obtained 
by substituting Eq. ( 1 , 4 _ 15 ) in (1.6-9) for q , i.e., 


Po 9 ^ » e) + Pq q + Q^l i + = 0 (1.8-1) 


where 


w = T^i) ^ + }(ko)i y 

d-^ -^-ij d^ ilok 


and represents the dissipative work done. Since, 


(0 Q^)|i = eji + ee'-ji 


( 1 . 8 - 2 ) 


(1.8-3) 


multiply Eq. (1.8-1) by 6 and then substituting the expression 
for eQ^,. from (1.8-3) 


pQ ee ^ (^ , G) + p^ q 0 - 01 . Q^+0 + (0Q^)i . = 0 


‘d "11 

(1.8-4) 


Since this is true for the whole volume or a part there of, 

integrating the expression over the volume v and then trans- 

* "1 

forming the volume integral of the teim (0Q ) j ^ to the sur- 
face by G-reen-Gauss theorem, Eq. (1.8-4) yields, 


/ ( Po 99 ’ 9) + P^ 0 q - 0j^Q^ + w^e) dv 

+ f n. eQ^ ds = 0 (1.8-5) 

s 

which is to be associated v/ith the relation on the boundary 
surface, 


h = 


n^Q 


i 


( 1 . 8 - 6 ) 
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Here s is the 1100116.3x7 surface enclosing the volume v and n^ 
is the unit outward normal at a point X. In Eq. ( 1 . 8 - 5 ), 9 
is temperature, is heat flux vector, q is the distributed 
energy source inside the body and ^ is free energy and in 
(1.8-6), h is heat flux. Constitutive equations are to be 
supplied for 4 and Bq. ( 1 . 8 - 5 ) forms the basis equation 
for thermal state in integrodifferential form. It will be 
kept as it is for later derivation. 

1,9 Constitutive Equations: 

The general physical laws in themselves do not suffice 
to determine the deformation of motion of a body subjected 
to given loading. Before a determinate problem can be for- 
mulated, it is usually necessary to specify the material of 
which the body is made. In the discipline of continuum 
mechanics such prescriptions are designated- as constitutive 
equations. The derivation of constitutive equations is a 
fundamental issue in the theoretical development of nonlinear 
mechanics. The basic derivations though based on much mathe- 
matical rigour, are not unambiguous regarding their fxmdamental 
11 12 20 

assumptions * ’ . Moreover, from practical point of view, 

the constitutive equations are not exact descriptions of even 
the gross physical properties of real materials. The best 
that can be said of any mathematical model of material beha- 
vior is that it provides a useful description of certain 
features of the behavior of some real materials. 
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The concepts of irreversible thermodynamics were 
applied to the study of linear materials v/ith memory by 
Biot*'^ and Ziegler^"^. These applications have been extended 
to include nonlinear materials undergoing large deformations 
by Talanis"'^, Ooleman and Gurtin"^^, and Coleman andNoll"*'^, 
Goloman and Guitin have pointed out that this approach to 
continuum thermodynamics is but one of several other appro- 
aches including those based on constitutive equations of 
differential type and on the axiom of fading memory. The 

differential type constitutive equations have been studied 

i ^ 21 

by Coleman and Mizel , Mandel and Brun and others. On 
the other hand, the axiom of fading memory was developed by 
Green and Rivlin^'^, Noll^^ and Coleman^^ and elaborated by 
Coleman and Mizel^^, Mizel and Wang^^ and others. The equi- 
valence of these two approaches, in certain special forms, 

2*7 2R 

has been studied by Coleman and Noll , Lubliner and 
oq 

Lianis 

Constitutive equations prescribed in the above lite- 
ratures are exceedingly complex and general. These equation 
are not usually suitable for use, in their full generality, 
in the solution of initial and boundary value problems. The 
purpose of such equations is to provide a framework within 
which observed behavior under a wide range of conditions can 
be correlated, and to serve as sources of simpler equations 
valid under more narrowly prescribed conditions. Within 
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such a framework various types of seemingly contradictory 

50 ■ 

special equations can he reconciled , as narrow range 
approximate descriptions. However, in this variant not much 
of stress will he laid on the mathematical basis for the 
derivation of the constitutive equations and on their res- 
trictions and the physical significance; on the other hand, 
fairly general forms of stress, strain and temperature laws 
will he prescribed which can he simplified for particular 
cases. However, various special forms of constitutive equa- 
tions have been listed in the reference. 

The t-hermodynamic stresses and dipolar stresses 

51(^0 )i governing equations (1.7-5) and (1.7-6) are of 

two parts; elastic or recoverable and dissipative or irre- 


coverable. The principle requirement is that they should 
not violate the laws of thermod3niamics and principle of 

''Mil 

material frame indifference. Recalling that T' 'is iden- 
tically zero, has been assumed to be nonexistant and 

that entropy can always be written as a function (or more 


generally as a functional) of the histories of 9, and 


E. and their time derivatives upto the present time, the 

1 J K. 

following constitutive equations can be prescribed: 


rj ( 0 i) =5} ^Pl) (E 


9, 9, tP) ( 


dU. , . 
3 1 
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(31)= T^pq) 


d-" ^\l’ \lm’ \l’ \lm’ ®’ tP)(gy^^ + g^^^- )/2 


6E ai 

■.|3 ^^3|i- 

( 1 . 9 - 2 ) 


j(ko)i^ g^(pq)r fE E A A o a fP^ 

r ^klm’ ^kl» -^liii’ ©» e» ^ i 


8 E 

Pn rqp 


8U. 


ijkj 


In relations (1,9-1) and (1.9-2) the argiment t^ denotes the 
dependence of the ; functions on the history of their arguments 
for t > t^ > •- 00 , where t is the present value of time. The 


arguments A's are Rivlin-Ericksen tensor 

these tensors are given hy, 

(p) r pm(a) (r-a) 

A. . = Z (^) (x Y I . V I . 


3 


Eor a general case, 


(r) r 

and A. = ' s 

cc=1 “ 


( 0 ) 

with Y I . = & . 

m 1 mi 


( ) Cr^ 
a'^ 

( 1 ) 


(r-a) (a) 

Y i . Y I 
mi nhk 


( 1 ) 


1. . = A,, and = A. 


( 1 . 9 - 3 ) 


In (1.9-3), the hraces are permutation symbols 




r! 


a: 


( r-a) ! 


(r) 


(r) 


and Y^j ^ "'^ml 3 k higher order velocity gradients i.e., 


(r) 

Y 


m 1 


E^(U I . ) 
^ ml 1 ^ 

Dt^ 


(r) 

““I \| jk 


“Xijik 


Dt- 


In constitutive equations (1.9-1) and ( 1,9-2) only the velo- 
city gradients upto one (r = 1 ) has been considered directly 
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ajid any higher order effect can he incorporated through the 
definition of the fimctional T(.) and its histories. Since, 
in the next chapter, an incremental approach will he adopted 
for Solution which v/ill actually move through the range of 
the parameter (history), no special consideration need to he 
attributed to these functionals (for an elegent proof see 
reference 31). 


For example consider the stress strain law of a visco- 

22 

elastic material formuilated on the basis - of n step history 
Since an arbitrary history can be approximated as closely as 
desired by a step history, the stress components which can 
be written in the form, 

^ oo 

(t) = I ^(Pl) (t) (1.9-4) 

n=1 ^ 


v^rill give the exact stress for any history, assuming reason- 
able smoothness of dependence of stress on the strain history. 
Pipkin and Rogers^^ have taken (t) in the form of multi- 

ple integrals which, in present notations for three dimensional 
case can be shovra to be, ■ 


T 


(pq) 


n 


(t) = m -t 


.... Ei.(t*), t-t* ). (1.9-5) 


where, 


6 ,' = 

k n ^ 


3R 


(pq) 


'n 


m 


rs 


— hs p « 


k 


k = 1,2, . . . , n 


(1.9-6) 
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S’or one step method, it takes the simplified familiar form. 




T 




— OO 


rs 


rs^“V 
(1.9-7) 


Again consider a non-viscous elastic plastic medi-um 

having no memory (a more general case for a medium with 

’ 5*5 

memory has been treated by Kluitenberg ) . In this case, it 
is possible to introduce the physical assumption that if 
plasticity phenomena occur, it will be governed by the yield 
function ^ , defined in cartesian coordinate system, by 

$ = ~ ^ (1.9-8) 

i,D = 1 

which is Von-Mises yield function. The associated Von-Mises 
flow law can be expressed as. 





(1.9-9) 


where b vanishes for ^ < 0. for such substance, the 

stress strain relation is given by 

^5.(13 ) 


ab* T 


(id) 


dt 


= a 


dt 


( 1 . 9 - 10 ) 


where a is a material constant and b vanishes for ^ < 0. 

This is Prandtl-Reuss equation. ■ 

In the equation for thermal state (1.8-5), the cons- 

* 

titutive unlcnowns are t and Q^. from the similarity of 
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expressions (1.7_3), (1.5-15) and (1.5-10), tbe fnnctional form 
of the Helmholtz free energy, can he assumed as, 


/V 

/s. 



which for a linear elastic material 


^i^k’ (1.9-11) 

, can be explicitly written 


as , 

p ^ = -L o^3kl 2. . E, ^ +0^3 E. . 9 + 1 U 9^ (1.9-12) 
^0 2 13 kl ~ 13 2 — 

where q material coefficients, for linearly 

elastic case, these constants may he assumed to he the function 
of 9. For nonlinear material they can he considered as func- 
tion of 6 as well as j • The coefficients and , 

have the S 3 nmiietry restrictions, 

Qijkl _ Q3ikl ^ Qijlk _ Qkli3 

and 9^3 = (1.9-15) 


As regard the heat flux vector Q^, a functional consti- 
tutive equation can he written in the form, 

. . 


Q 


1 _ 


^®kl’ \lm’ 


®!k’ 


(1.9-14) 


Coleman et al. have proved through the restriction of second 
law of thermodynamics that for a rigid 'heat conductor, the 
classical Fourier’s lav/ of heat conduction can he obtained 
as a first order approximation having the form, 
qi = Q| 

K^3 = e:^3(0) 


(1.9-15) 



I»l. T. KANPIM 

C£!Vi<Ui imAHY 


I jiiS87s^ I 

1,10 Summary : 

1.10.1 Assumptions: 

The following assumptions have “been made in the preced- 
ing sections for deriving the general governing equations for 
continuum. 

1 . The material will ohey conservation of mass and first 
and second laws of thermodynamics. 

2. Hypothesis of m.aterial frame indifference is applica- 
ble. 

5. Concept of force defined through work rate principle 

is valid. 

4. Displacement derivatives upto second order will he 
active for the internal energy expression. 

5. Supply of nonmechanical energy is only due to thermal 
effect. 

6. Multipolar kinetic energy and heat flux vector has 
been neglected. 

7. Antisymmetric parts of the dipolar stresses, Ei.e., 
t! ^3 ii nonexistent. 

1.10.2 Governing Equations and Constitutive Relations: 

The governing equations for the continuum can be 
rewritten from Eqs. (1.7-20), (1.7-21) and (1.8-5) for 
mechanica.1 effect as, 

f (T^^i . + p _ P yi) Y, dv + / (P^-n.T3^)Vi ds = 0 

V ^ ^ ^ ^ ^ s 3 

( 1 . 10 - 1 ) 
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and /(rpDi_rp(oi) + ^(kD)i + p pji) y 

V I k; 0 1 1 0 

+ / T(kd)i) ds = 0 (1.10-2) 

s 

and for thermal effect as, 

/ (Pp 0 e ^ (ij^, 0) + Pp 0q - Gj^ + w^6) dv 

+ / n. 0 ds = 0 (1.10-5) 

s ^ 

where w, = A.. + ^(kj)! ^ (1.10-4) 

dd r^d i-lDh 

Here w is any arbitrary volume ho’unded hy the surface s having 
an outward normal n^^ at a point X. P^, and P^^ 

in (1.10-1) and (1.10-2) are respectively monopolar stresses, 
mono and dipolar body and surface forces. T^3i) ip(h3)i 
are thermodynamic tensions and is the velocity. In (1.10-3), 
0, q and are respectively temperature, distributed heat 
source . and heat flux vector. The symbols and p^ stand 
for dissipative ?7ork rate and density of mass, while i> denotes 
free energy, ill the quantities referred here are in terms 
of reference configuration. The heat equation (1.10-3) must 
be complemented with the relation, 

h = n^ (1.10-5) 

on the surface s, where h is the heal flux over the surface. 

Specification of mechanical and thermal properties of 
matter through constitutive relations for mono and dipolar 
stresses, free energy and heat flux vector may be prescribed 
through any general functionals provided they satisfy frame 
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indifference principle, laws of 'tbermodynamics specially 
second law restrictions and any other physical character- 
istics imparted on the material properties. Also, as per 
asstunption, they are not supposed to contain derivatives 
of displacements of order more than two. In restricted 
sense, they are prescribed hy Eqs. (1.9-1) to (1.9-4) and 
are given as follows: 
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The constitutive relation for free energy’ll has been 
set as functional, 

^\m’ ®lmn’ -^m’ ^"imn’ ®’ ®’ (1.10-7) 

and heat flux vector, in a simpler form, as 


(E^j^, E^^, 0.^, tP) 


( 1 . 10 - 8 ) 


In the above equations, E .^^ , and may be chosen 

as any strain measures, but conforming to tradition, only 
classical measures have been quoted here. 
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With slight modification, and can he xvritten in the 

form, 
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The symbol t^ has been introduced to show the dependence of 
the functional on the past history of the argument s. Hence, 
if t is the present time, then t > t^ > - . 


In the next chapter, finite element method has been 
discussed and the above set equations have been transformed 
for finite element use. 
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CHAPTER 2 


DISCEETIZITION MD inffliRISATlOH 

2,1 Introduction: 

The abysmal mathematical complexities encountered 
for the exact solution of differential equations for conti- 
nuum mechanics are virtually insurmountable except for very 
few idealized cases . Even the classical approximate pro- 
cesses such as as 3 nmptotic expansions and weighted residual 

2 

techniques can be applied to relatively simpler problems , 

In the past few decades, due to the introduction of high 
speed Computers, some classical approximate methods such as 
Ritz method, finite difference technique, weighted residual 
method have been revived and extended. Development of new 
concepts like invariant imbedding, finite element and other 
numerical procedures have also been introduced. 

Although, finite difference method has been mathema- 
tically studied rather well, there are inherent difficulties 
associated with its application in continuum mechanics prob- 
lems. The inability to express some of the boimdary condi- 
tions and the unaccommo dative nature of the method for irre- 
gular changes of geometry and material properties are some 
of these difficulties. Eor the method of invariant imbedd- 
ing^’^, which is an extension of Rewton-Raphson technique 
for functional analysis, the disadvantages are difficult 
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formulation and instead of one original problem a family of 
problems must be solved resulting in more computational effort* 
imong the approximate techniques, the finite element method 
can be considered as one of the most powerful and versatile 
discretization presently available for the numerical matrix 
solution of complex structural problems using digital com- 
puters. 

In conventional matrix structural analysis, the system 
is usually composed of a finite number of structural elements 
connected only at a finite number oi nodes, Ihe stiffness 
(in some general sense) properties of the elemental subdivi- 
sions can be evaluated exactly within the scope of the adopted 
theory. Depending on the variables used, in the definition of 
the stiffness matrix, the method can be classified under three 
heads - (a) displacement method, (b) force method and (c) mixed 
method. 

In finite element method, the domain of the continuum 
is discretized into a number of disjoint subdomains called 
elements represented by a set of points which are called nodes, 
Ihe elements may not only be connected at the nodes, but also 
along the interelement boundaries, Then the unknown field 
variables in the integrand of the functional integral repre- 
senting the equation of motion of the system, is approximated 
by a set of assumed functions which are expressed in terms of 
the field variables at the nodes by suitable interpolation 
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formulas. The stiffness (in some general sense) properties 
of such elements can now he defined approximately and can 
he determined in many alternative ways. By far the most 
popular method for obtaining finite element equations is 
based on variational principles with respect to minimization 
of potential or complementary energy or Reissner’s energy formu- 
lation. However, there is no potential harrier of applying 
other approximate procedures such as Galerkin’s technique or 
any other principles in mathematical physics. 

In the problems of continuum mechanics, the finite 
element method is easily adaptable to matrix formulation which 
can be readily analysed by digital computers. The method is 
capable of approximating quite irregular boundary shapes, com- 
plicated boundary conditions and arbitrary, variation of phy- 
sical and geometrical parameters in its domain of applicabi- 
lity, Here, the algebraic equations resulting from finite ele- 
ment discretization will be, in general, nonlinear and simul- 
taneous, In the absence of any direct procedure for solution, 
some iterative and/or an incremental technique has to be adopted. 
The combination of the finite element method with incremental 
forward integration procedure on the field variables is quite 
suitable for computer use and general problem solution. 

In this chapter, after giving a brief historical review 
on finite element in section 2 and types of models used for 
analysis in section 3, different modelling criteria including 
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tlie question of conirergeucy is discussed in section 4, The 
crux of the finite element analysis lies in choosing the dis- 
tribution of the field variables. In section 5, the field 
variables have been chosen in a general form without any 
restriction on the Interpolation law and the geometry of the 
element. In sections. 6 and 7, these distributions have been 
substituted in the equation derived in Chapter 1 to get cor- 
responding equations for finite elements. In section 8, 
parametric differentiation scheme has been introduced for 
stepwise solution. 

2.2 Historical Review: 

The development of discrete method of structural mecha- 
nics had its beginning in the early 1950*s together with the 
advent of high speed computers. The concept of finite element 
was first introduced by Cour^t"^. He solved in an approxi- 
mate manner the classical St. Tenant torsion problem by assimi- 
ing linear distribution of warping function over each triangular 

doittains and then by applying variational principle leading to 

6 

minimization of potential energy. Prager and Synge and 

7 

Synge' further elaborated the technique by a geometric repre- 
sentation in function space which is often called hypercircle 

Q 

method, Argyris * text for the matrix formulation of the 
transformation theory of structures is an important milestone 
in this development work, facilitating a better understanding 
the structural behaviour and providing a practical and powerful 
way of employing automatic computation. 
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This method was first applied to plane problem hy 
Turner and his associates^ using triangular and rectangular 
element model. Very soon after its first application, it 
"became evident that the main feature of the displacement 
consistent finite element is the assumed kinematic field 

distribution and the method extended to three dimensional 

10 11 12 
analysis , plates , shells of revolution , axisymmetric 

bodies shallow shell , general shell structure eccen- 

1 6 

trically stiffened cylindrical shell with considerable 

amount of success in each case. During the past few years, 

study of the mathematical foundations of the method and 

1 7 

different types of variational principles applied for the 
■ generation of more general, and varied types of finite element 
models as well as its applications » have greatly 
extended its fields and clarified the basic requirements for 
its effective formulation. 

In early stages of development, a theory establishing 

necessary and sufficient conditions to guarantee convergence 

to the true solution was lacking. Such requirements were 

2 "^ 

first considered by Melosh . The development of the stiff- 
ness relations may be shown to be equivalent to piecewise 
Rayleigh-Ritz procedure applied to variational principles, 
provided continuity requirements at the interfaces of the 
elements are satisfied. The fundamental considerations for 
the successful generation of stiffness rest on the choice of 
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distri'b'ution of unlmo'wn field variables having the require- 
ments of (a) isotropy, (b) completeness and (c) continuity 
conditions which have been discussed in references^^*^^’^'^. 

The failure of early attempts^^ for the satisfactory construction 
of triangular elements for plate bending analysis may be attri- 
buted due to ignorance of these criteria. 

With a firmly established mathematical basis for the 
generation of finite element stiffness matrices, a systematic 
development for refinements is possible which may be observed 
in the lst^^ and 2nd^"' conference on Matrix Methods in Struc- 
tural Analysis, The application of such elements to geometri- 
cally and physically nonlinear problems which were previously 

restricted to framed structures, have been possible with consi- 

53 54 

derable confidence and certainty Time dependent pheno- 

mena^^ and plastic analysis^'^»^®»^^ have also been done 
employing this technique, 

finally, it should be emphasised that finite element 
technique is by no means restricted for the solution of struc- 
tural problems. It can be, as well, applied to any field prob- 
lems with considerable ease, Zienhiewicz and Cheung^^ have 
applied this procedure to solve second order problems such as 
stationary heat flow and torsion, while application to coupled 

thermoelasticity, nonstationary heat conduction, fluid mecha- 

51 

nics problems may be found in reference , 


It may be noted 
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here that hhe references cited above are by no means exhaus- 
tive which is mainly because of the ever expanding litera- 
ture on this subject, 

2,3 Types of Model: 

Depending upon the characteristics of a finite ele- 
ment model, two principal classes can be defined: (a) dis- 
placement model which inherits the chosen displacement field 
and (b) eq_uilibrium model which starts with the assumption 
of an ‘eq^uilibrated stress field. Clough, Argyris, Zien- 
kiewicz and their associates are among the' notable contri- 
butors on displacement model, Veubeke^"^ and Morley^^ have 
applied equilibrium model to plate bending analysis. The 
displacement and equilibrium models may be considered as 
subsets of ’mixed model’ where both displacements and stresses 
may be selected independently. Herrmann^^, Prato"^^, Dunham 
and Pister have developed mixed models for various indivi- 
dual problems employing Reissner’s -variational theory. Again 
depending upon the method of solution there may be three 
classes, (a) displacement method, (b) force method and 
(c) mixed method. However, such classifications are not 
•unambiguous^^, Por example, there may be three variations 
of finite elements which finally employ nodal displacements 
as unknowns. The first variation starts by choosing displa- 
cement distribution in the individual elements and equations 
are formulated by application of minimum potential energy 
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principle, This is a displacement model. The second varia- 
tion uses an assumed equilibrated stress field and derives 
the basic equations from the principle of minim-um complement- 
ary energy. The stress field is subsequently eliminated in 
terms of conjugate generalised nodal displacements. This 
may be designated as equilibrium model. The third variation, 
a hybrid model as has been designated by Pian, is based on 
modified complementary energy principle, for which compati- 
ble displacement functions are assumed along the interele- 
ment bounda,ries, in addition to the assumed equilibrated 
stress field in each element. The above three variations 
may be classed as displacement method. Similar subdivisions 
can also be made in the force method. Hence, a dual hybrid 
model can also be formulated ass-uming equilibrated tractions 
along the interelement boundaries and continuous displacement 
field in each element. More exhaustive discussions, examples 
and references are given by Pian and Tong in addition to 
the associated variational formulations in each case. 

2,4 Modelling Criteria: 

2.4.1 Efficiency: 

in attempt to assess the usefulness of a finite 

element model must be based on some efficiency criteria, 

24 

Poliowing are the three criteria that seem to be generally 
acceptable. 
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1 . Convergence : 

As tlie sizes of elements in a model approach zero, the 
discretization error in computed field variables, as well as, 
in geometric simplification, if any, should approach zero. The 
convergency should be monotonic, 

2» Rate of Convergence: 

The efficiency of a model should be proportional to 
its rate of convergence, 

5. Computation Efficiency: 

The efficiency of a model should be proportional to 
the ratio of accuracy to computational effort. 

Depending upon the above characteristics, there should 
be some way of expressing the efficiency of a given model from 
which usability of different models can be compared, 

2,4,2 Mathematical Criteria: 

It has been already mentioned that the distribution 
functions of unloiown field variables should satisfy the xe- 
quifements of: 

1, Invariance, 

2, Completeness and 

3, Continuity. 

In early stages of development these requirements were 
not systematically considered which has resulted in many 
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Similar convergency criteria were discussed 

px pc 

on heuristic "basis "by Melosh , and Irons Johnson and 

Mclay"^^, Tong and Pian^"^, Walz, Pulton and Gymis^"' have also 

investigated convergence in connection with the particular 

cases of finite element models. Oliveira’s works concerns with 

27 

the mathematical "basis and in particular, completeness and 

2B 

convergence criteria . He has derived the sufficient condi- 
tions for convergence of a general case of displacement ana- 
lysis with no particular reference to any definite model. 
Although, the analysis is restricted to displacement model, 
it appears that similar conditions can also he applied to 
any other model. 

In classical mathematical formulation considering va- 
riational or integral equation approach, it is usually required 
that the assumed functions should possess derivatives which 
are continuous upto highest order occurring in the correspond- 
ing functionals except possihily at isolated singularities 
like points, curves or surfaces which necessitate special 
attentions. . In finite element analysts, the distribution 
functions satisfy the strict continuity requirement within 
each elem.ent, but on the interelement boundaries the admissi- 
bility conditions are relaxed in such a way that the functional 

19 

will he defined for the chosen distributions 

1 , Invariance : The invariance condition will guarantee that 
any rigid body parallel shift or rigid body rotation of the 
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coordinate system will not change the computed field variahles 
except that they are now the components in the new coordinate 
system. This shows that the chosen distribution should satis- 
fy the invariance requirements given in sec, 3 of Chapter 1. 

In the case of polynomial distribution function, it is necess- 
ary to use complete polsraomial upto a certain order, since it 
is isotropic and obeys the frame indifference principle. In 
the Case of natural coordinate system (which is one having 
the element interfaces as coordinate surfaces, e,g. , area 
coordinates), the invariance requirement will be automatically 
satisfied if the distributions are represehted by isotropic 

nc 

functions of the natural coordinates , 

2, Completeness: 

The satisfaction of completeness requirement will 
ensure the convergence of the field variables to the exact 
values as the element sizes approach zero, provided the con- 
tinuity conditions are not violated. However, convergence to 
the true solution may still be achieved even if the continuity 
conditions are not satisfied, whenever the body fo3XJe density 

corresponding to the approximate solutions remains bounded 

28 

within the element, no matter their size , 

Completeness of a sequence of families = 1,2,...) 

with respect to a given set { UjJ is only meaningful provided 
it is possible to establish a predefined norm relating 
and each member Uj, of the setiUj,}, Then the sequence 
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will be designated as complete, if for a specified 6 > 0, 
it is possible to obtain an integer N such that there exists a 
member v/hich satisfies, 

II “ij. % II < e 

for any IT > N and for any member U-g of { Ug) where j | . j ] 

is the predefined norm measure. 

Unfortunately, it is not always easy, in general, to 
prove directly the completeness of a distribution function. 

But as has been mentioned in'^ ’ , for a, polynomial distributi„ 

completeness will be . achieved provided, 

(a) the distribution function has a number of arbitrary 
parameters equal to the number of independent unit modes 
corresponding to the element; 

(b) the terms of degree higher than first can be varied 
independently with respect to the constant term and the coeffi 
cients -which affect linear terms; 

(c) the constant terms and the coefficients which affect 
the linear terms are completely arbitrary. 

When using incomplete polynomial expansions in natural 
coordinates, it is not immediately evident whether the expan- 
sion is complete or not. A simple hut indirect way has been 

^ 26 
discussed by Belippa . 
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3. Continuity: 

The requisite interelement continuity of chosen field 
Yariahles, does not affect the convergence of the solution, 
hut it affects the rate of convergence. Lack of continuity 
may cause non-monotonic convergence of the solution. The 
method of obtaining the necessary continuity conditions to 
be satisfied! for a particular model have been discussed by 
Pian"*^. 

2.4.3 Geometrical Criteria: 

In finite element -method, discretization is not only 
done on chosen field variables but also on the continuum 
itself, i.e,, a curved member may be approximated as an assem- 
blage of small triangular flab elements. It has been shown 
by Walz et al that convergence may not be achieved if a 
curved structure is discretized by series of straight elements, 
but the convergence may occur when the element are curved. 

Apart from the error resulted from the mathematical 
or geometrical discretization, there may be round-off error 
associated with the accuracy with which the numbers are 
manipulated in a computer. Por solution of a problem having 
a large number of elements, the error resulted from this 
source may also be critical. However, round-off error is 
not the characteristic of a finite element model, though 
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it indirectly effects the solution depending on tbe order of 
interpolation formula used to describe the field variables. 

In the subsequent sections, finite element equations 
will be obtained for a general case without restricting to 
pay particular model. 

2.5 Finite Element Derivation: 

For the construction of finite element model of the 
basic field equations (1.10-1), (1.10-2) and (1,10-3) of 
the last chapter having independent field variable as displace- 
ments U^, stresses and temperature 0, consider the con- 
tinuum to be represented by R number of elements with a typi- 
cal element as m and typical node as p and the number of 
nodes surrounding the node p as M. Let .a^ , Ugj aJid 
be the number of degrees of freedom per element for the inter- 
polation of displacements, gradients of displacements, stresses 
and temperatures respectively. The interpolating functions of 
the field variables for a typical element may be approximated 
as follows ; 

1. Displacement distribution: 

The displacement distribution over a typical element 
m may be prescribed in. terms of generalized nodal defomna- 
tions as, 

(X^, t) = (2.5-1) 

(X i j ♦ • * 5 cx.| 
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where k takes the val-ues of one, two and three representing 
the three space variables. The generalised deformations need 
not he nodal displacements only, they may imply gradients of 
displacements as well. The contravariant components of dis- 
placements can be represented as, 

U3(X^, t) = D.^(X^) U^(t) (2.5-1a) 

where is the contravariant metric tensors and are 

physical components of the generalised deformations at the 
nodes. 

2, Displacement Gradients: 

Dor convenience, the gradients of displacements can be 
approximated within element m by the corresponding generalized 
nodal gradients Rppf(t) by, 

hi 3 = hjp 

In equation (1.10-2), the constitutive expressions for T^^^^ 
and involve gradients of the displacements. These 

gradients are evaluated by differentiating Eq. (2.5-1). Thus 
they can be obtained in terms of the nodal displacements only. 
This means that the derivatives of the nodal displacements are 
not involved in the constitutive equations. However, (2,5-2) 
may be used for the displacement gradients in Eq. (1,10-2) in 
the last chapter. Hence Ilgj;^(t) is assumed to be independent 
variable and not function of the nodal displacements. 
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3. Stress Distribution: 

The distribution of stresses can be similarly assumed 

to be, 

1"^=^ (X^, t) = (X^) Tpi^xCt) (2.5-3) 

P *1 y • * • 5 (X^ 

where fpjj(t) are generalized nodal stresses and are 

corresponding interpolation functions. 

4. Temperature distribution: 

In exactly similar way, the temperature distribution 
can be approximated in terms of generalized nodal temperatures 

0(X^S t) = (f") 0^jj(t) (2.5-4) 

Y "I,..*, 

The field variables T^^j and need not be con- 

sidered as tensor components, but the summa-tion convention for 
the repeated indices has been used in the above expressions 
and will be retained throughout, unless otherwise stated. 


2.6 Element Equations: 

Introducing distributions (2.5-1) to (2.5-4) and 
applying summation extending over all the elements, Eq. (1.10-1) 
reduces to 


R 

S 

m=1 




T„,t - 02 , 

pH aa ’ 


a’H 


11 ) = 


0 


( 2 . 6 - 1 ) 
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where a, a’ = 1 , . . . , ; P = 1 , . . . , 

°hp = { h. "j ha®f "" 

in m 

02 , = / p ID. , dw (2,6-2) 

aa* :l la a' 

m 

L1 = / P ^'^dv + r D. ds 

a t '^o la “ la 

V b 

ra E 

The limit correeponds to the volume occupied hy the ele- 
ment n hounded hy the surface s^. In Eq. (2.6-1), will 
he common only to M surrounding elements among the total number 

R in the continuum. Again, since Eq. (2.6-1) remains true for 

« 

any arbitrary (j^O) , it implies 




m=1 




aa' 


U 


a ’IT 


II ) 


0 


(2.6-5) 


where the summation E* extends over the surrounding E elements. 
Sq. (2.6-3) provides the governing equation corresponding to 
(1.10-1) in finite element model. Similarly, substitution of 
distribution (2.5-3) and (2.5-4) in (1.10-2) and (1.10-3) of 
Chapter 1 , will result 


M 

° (2.6-4) 

in=1 

M 

and 2’ (04 , 6.^,^ - 14^ + 15^ + ^6^ ) = 0 (2.6-5) 

m=1 YY. ^ ^ ^ ^ . 
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16 = / S “a i’' 

T Y 

n 

It may "be noted that (2,6-3) provides three eq'ua-tions 
for each node, corresponding to displa-cement interpolation. 
The ntmiher of eciuations from (2.6-4) is six or nine pet node 
of gradient interpolation, for the monopolar or dipolaf cases 
respectively. Only one equation per node (for temperature 
interpolation) arises from (2.6-5). 

In each element, the total number of equations from 
(2.6-3), (2.6-4) and (2.6-5) is + Ug + a^, where as the 
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total n-uinber of imlcno\sms are + a-? + a. . Hence, for deter- 
mlnancy of Solution, a 2 Has to He equal to a^, 

The set of-Eqs. (2.6-3) to (2.6-5) is nonlinear coupled 
simultaneous differential equations with respect to time t. 
Einite element transform.a.tion hs,s reduced the continuous fimc- 
tion of field variables , ’with respect to space coordinates, to 
a finite number of field unlsnovms rendering the system to be 
space independent. This reduces- the complexity of the non- 
line a,r equations and become more amenable for solution by 
computer operation at the expense of increased number of 
governing equations. 


2,7 Expressions Eor Strains: 

The above equations in the foregoing section are not 
complete without explicit functional relations of E. . and 

i- J 

E . ... It has been alrea.dy mentioned that various expressions 

1 JX 

may be chosen for straans, their principa.1 requirements being 
that they should be frame indifferent and may be treated as 
arguments of the internal energy expression, reflecting the 
physical behaviours of the material-. Eor the saEe of brevity, 
the strain measures given in (1,10-9) to (1.10-12) are pre- 
sented here. Eor future usage, the following quantities have 
been expressed in teims of displacement distribution. 


1. Infinitesimal strain: 


B. . 

ij 


-m.iP/S =4-0 


IN 


a 


ON 


(id) ^'cdT 


(2.7-1) 
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2. Infinitesimel rotation; 


^13 




E. . = E. . + (E . + iR .) (i . + 


3. Classical Strain: 

13 ID 


where 


mi 


’mi' 'no no 


= 

1 

“2 

(0°^(io) ^ ^Tio) ^oll 

(2.7-3) 

^7 13) 

z= 

D . 1 . -S' E . 1 . 
ia|0 Oaji 

(2.7-4) 

^7 13 1 

= 

E. , . - E. 1 . 
la] 0 Docii 

(2.7-5) 

^7^3 ) 


nmn 

V(^7ml) ^7i>3)"^7ml) 

^ Uo ! 


+ I l) (2.7-6) 


(no) 


mi 


no 


and G™' is the metric tensor in undefomed state. 


4. Dinolar Strain: 


E. 

lOk 


( r, 


i(ok) 


+ ^ 


( aa’ ) 


U ,^) u 


iCok) a'N''^ oN 


where 

^(Slc) 


^iaj ok 

G™ D E 

maji na’jok 


(2.7-7) 

(2.7-8) 

(2.7-9) 


It is evident from (2,7-6) and (2.7-9) that if some nonlinear 
terms in E^^ and E^^^ are to he neglected (i.e., assuming them 
to he small in comparison to others due to some physical 
reasoning), they may he easily incorporated in the expressions 
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for smd 'the simplified expressions, care 

AR 

should be taken that they represent integrable functions^ . 


5. Rivlin-Ericksen Tensors: 



(2.7-10) 


A 


ijk 


0 “ 


i(Dk) 



(2.7-11) 


In all the expressions from (2.7-1) to (2.7-11), a and a* 
vary from 1 to . 

Although, the finite element equations are simpler 
than the original set, they are not so simple as to be amen- 
able to straight forward solution, unless they represent a 
linea,r case. In the absence of any direct and definite proce- 
dure, the method of parametric differentiation technique has 
been employed to obtain a step-by-step solution. 

2,8 Parametric Pifferentiation: 

2.8.1 Introduction: 

The solution field of a system may be regarded as a 
phase space, being functionally related with a forcing para- 
meter causing the motion of the system from one phase point 
to another. The philosophy behind the technique of parametric 
differentiation Is to evaluate at a certain phase point, the 
tangential stiffness of the system for its motion due to a small 
perturbation of the forcing parameter. According to this view 
point, the method may be described as a generalized form of 
either infinitesimal perturbation or incremental approach. 
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The technique of parametric differentiation was first 
discussed hy Davidenko^^. At a much later date, the charac- 
teristics of this method, was discussed hy Yalcolev^*^, Appli- 
cation of this method for the solution of different problems 
may be found from the works of havidenko^"* , Setlur and Kapoor^^, 
Jischke and Baron^^, Rubbert and landahl^^. The method essen- 
tially consists of choosing a parameter which may be assumed 
as an independent variable X for all the field variables and 
for which the solution is known at some limiting initial value, 
say at X ~ ^ o * nonlinear equations with the boundary con- 

ditions, if any, are then differentiated with respect to X to 
yield a system of linear equations with variable coefficients in 
terms of differentials of the field variables and may be solved 
employing any standard numerical integration technique. The 
solution of the problem can now be obtained by further inte- 
gration of the differentials for an assumed step size at X^, 

The Solution can thus be marched out from X^ to any other 
value of X . The method can be best demonstrated by taking a 
simple example. Say it is required to solve the nonlinear 
differential equation, 

CU + 1 = 0 (A) 

with boundary condition 

U = U, at t = t^ (B) 

X o 

o 

♦ • 


where 


0 = C(B, TT), B = I,(t), TJ = Tj(t) and tJ = dU/dt 
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Choose the forcing parameter X such that 

L = L(t, X ) where L(t) = h (t, X^) 
and U = UCt, X ) where U(t) = iJ (t, X^) 

Consequently, it is possible to construct 

c = C(ij(t, X ), Ti (t, X )) 

where, 

C(t) = U (u(t, x^), n (t, x^)) 

Hence the new differential equation can be written as, 

C . ^ + 1 = 0 (0) 

with boundary condition, 

U (t X ) = U , (D) 

O Oq 

The function U has to be chosen in such a way that it satis- 
fies (C) and (D) at X = X^. 'Sow differentiation of (C) and 

(D) with respect to X yields, 


3U 


+ ( ^ + C ) H' 
au 


0 


(E) 



The modified equation (E) and (E) is linear with respect to 
U and with variable coefficients. It can be solved, at least 
numericaJly, and the solution will be U*(t, x^). From this 
U(t, Xq + AX) can be easily obtained by the integration, 



83 


U(t, +AA ) = U(t, A^) + 

1 +AA ^ 

/ ° U (-fc, A^) dA (G) 

In similar. way, the solution can he marched upto A = A ^ 
when U will he the same as TJ which is the solution of (A) and 
(B). The advantages of this method are that it allows the 
solution to he obtained throughout the range of A and no, ite- 
rative technique is needed. 

Bor efficient execution of the quadrature in (G) the 
following two forward- integration methods are recommended: 

1, Improved Euler-Cauchy method. 

2, Runge-Kutta-Gill method. 

Both the schemes are easily available in any text of numerical 
analysis. 

2.8.2 Application to Continu'um Equations: 

Eor continuum mechanics problems, in general, there 
may be two t 3 ''pes of perturbations in the forcing functions, 
firstly, at a certain instant of tim.e, there may be jump in 
the loading, the time taken for this jump being negligibly 
small. Secondly, if the evolution of the phase points is 
considered as a smooth path with respect to time, the pertur- 
bations may be applied corresponding to small increment of 
time. That is, in the second case perturbation parameter may 
be chosen as time, where as, in the first case it may be 
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chosen as incremental loading. However, for the sake of 
generality, the parameter will he designated as A and all 
the field variables after finite element transformations, 
i.e., Ipjj- and Byjj including the body and surface forces 

will be assumed to be implicit function of A , Hence, the 
tangential stiffness will be related through equations which 
?/ill be obtained by differentiating Eqs. (2.6-3) to (2.6-5) 
with respect to A . It will be further assumed that all 
field variables are known at the starting point, A = A^. 

The object is to march forward from A^ upto a point defined 
as A = A^, which is also close to A^. By successive repe- 

titions, the solution can move to any desired value of ^ corres 
ponding to time or loading whatever may be the case. Hence, 
performing necessary differentiations on ( 2 . 6 - 3 ) to ( 2 . 6 - 5 ), 
the equation can be written in the form, 


M 

m=1 

M 

(C3pp, - Mp + 13p ) = 0 (2.8-2) 

m=1 

M 

y "X* 

and z' (04 , 9 ^ 6y«n - '^^y * ) = 0 

( 2 . 8 - 3 ) 

where a, a’ = 1 , . . . , a.| ; p, p' = 1, . . - , 

and y y ^ 
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In Eqs. (2.8-1) to (2.8-3), the asterisked quantities are the 
derivatives with respect to X of the corresponding unaster- 


isked quantities, t' 


U and 0.^ are respectively the 


cients, 




derivatives of nodal stresses, displacements and temperatures. 

These are unkno'vvns at this stage. LI and L3 n a 3 ?e incre- 

° a p 

mental mechanical loads and L5* is increment of heat due to 

Y 

increase of q. These are known quantities. The other coeffi- 

L2 p, L4^ and L6^ have to he expressed 

interms of the derivatives of the independent variables, 

* 

0 and their time derivatives. These quantities are very 
complicated and since solution of any specific problem with 
this generality is out of scope of this work, no attempt will 
be made to obtain the explicit expressions of the above quan- 
tities. However, it is not very difficult to visualize that 
the expressions can be obtained and after substitution, 

Eqs. (2.8-1) to (2.8-3) will be linear differential equations 


with respect to time and in terms of T 


* 


* 


pi’ ah 0 YN 
can be solved numerically. Once these quantities are obtained, 


the field variables 


, i.e. , T and 3 't X X 


'cdJ 

be ea,sily evaluated through the quadratures 

X 


1 


can 


U 


T 


oH 


pi 


IJ 


ccI 


X= X 


I 


0 


= T 


and 


^Yl 


pi 

0Y1 


X= X 


{ 




^ pi 


(2.8-4) 


X=: X 


o 


/ ' dx 

Xo ^ 
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The advantage here is that the dormant nonlinearity 
of the problem (since it is not (juite apparent in (2.8-1) to 
(2.8-3)) is confined only to first order equations and has 
formally been isolated in the quadratures (2.8-4) only. 

In the next chapter, application to a specific case 
concerning plane problems will be dealt with simplified 
assumptions on constitutive relations and explicit expressions 
of a3.1 the quantities will be obtained. This will further 
clearify the technique. 



87 


2.9 References: 


1. ^£ES, W. P. , Nonlinear Partial fPifferential Equations 
in Engineering, Academic Press', 1965. 

2. MffiHLIN, S., Approximate Methods for Solution of 
Differential and Integral Equations, iinierican Elsevier 
Publishing Co. , Inc. 1967. 

3. lEE, B.S. , Quasilinearization and Invariant Imbedding, 
Academic Press, 1968. 

4. BEIIMM, R. , KAliffiA, R. , and WING,G.M. , Invariant 
Imbedding and the Reduction of Two-Point Boundary- Value 
Problem to Initial Value Problems. Proc. Nat’l Acad, 
Sci., U.S. 46 (i960), 1646. 

5. COURANT, R. Variational Methods for the Solution of 
Problems of Equilibrium and Vibrations. Bull. Am. Math. 
Soc. 49, 1-23, 1943. 

6. ■ PRiGER, W. and SYNGE, J.I., Approximations in Elasti- 

city Based on the Concept of Function Space, Qrt. Appl. 
Math. 5, 241-269, 1947. 

7. SYI^GE, J.I, , Hypercircle in Mathematical Physics. 
Cambridge, Univ. Press. 1957. 

8. ARGYRIS, J.H. , 'Ener^ Theorems and Structural Analy- 
sis, Butterworths Sci. Pubs. , London, I960. 

9. TURNER, M.J., CLOUGH, R. W. , TLIETIU, H.C. and T0PP,L. J. , 
Stiffness and Deflection Analysis of Complex Structures, 
J. Aero. Sci. 23, 805-823, 1956. 

10. ARGYRIS, J.H. , FRIED, I,, and SCHARPF, D.W. , The 
Hermes 8 Element for the Matrix Displacement Method. 
Aero. J. 72, 691, 613-617, 1968. 

11. SLEYPER, H.A, , Development of Explicit Stiffness and 
Mass Matrices for a Triangular Plate Element. Int. 

J. Solids. Struct. 5, 3, 1969. 

12. ADELIMN, H.M. , CATHERIIIES , D.S. and WALTON, Jr. ,W.C. , 

A Geometrically Exact Finite Element for Thin Shells 
of Revolution, AIAA 7th Aerospace Sci. Meeting, N.Y. , 
Jan. 1969. 



89 


26. EELIPPA, C.A. , Refined Finite Element Analysis of 
linear and Non-Linear Two-Dimensional Structures, 

Ph.D. Thesis, Univ. of Calif. , Berkeley, 1966. 

27. OlITEIEjI, E.R, de A., Completeness and Convergence 
in the Finite Element Method. 2nd Conf. of Matrix 
Methods in Struct. Mech,, Wright-Patterson AFB, Ohio, 
Oct. 1968. 

28. OLITEIRjI, E.R. de A., Theoretical Foundations of the 
Finite Element Methods. Int. J. Solid and Struct. 4, 

10, 929-952, 1968. 

29. TOCHER, J.l. , Analysis of Plate Bending Using Tri- 
angular Elements, Ph.D, Thesis, Univ. of Calif. , 
Berkeley, 1963. 

30. 1st Conf. on Matrix Methods in Structural Mechanics, 
Wright Patterson Air Force Base, Ohio, 1965, AFFDL- 
TR-66-80. 

31. 2nd Conf. on Matrix Methods in Structural Mechanics, 
Wright Patterson Air Force Base, Ohio, 1968. 

32. ODEN, J.T., Finite Plane Strain of Incompressihle 
Elastic Solids by the Finite Element Method, Aero, 

Qrt. 19, 3, 254-264, 1968, 

33. WALEER, A.C., A Nonlinear Finite Element Analysis of 
Shallow Circular Arches. Int. J. of Solids and Struct. 
5,^ 2, 1969. 

34. lEMPNER, 0., Finite Elements, Finite Rotations and 
Small Strains of Flexible Shells. Int. J. of Solids 
and Struct. 5, 2, 1969. 

35. ODEN, J.T. , A Generalization of Finite Element Concept 
and its Application to a Glass of Problems in Nonlinear 
Tiscoelasticity. Div. in Theo. Appl. Mech. (Proc. in 
4th Sectam), Perg. Press. 

36. FRIED, I., Finite Element Analysis of Time Dependent 
Phenomena, AIM 7, 6, 1170-1173, 1969. 

37. RICHARD, R.M. and BIACKIOGK, J.R. , Finite Element 
Analysis of Inelastic Structures, AIAA 7, 3, 432-439, 
1969. 



90 


38. IS/iZSON, G, , Discrete Element Plastic Inalysis of 
Structures in a State of Modified Plane Strain Hill 
7, 3, 545-547, 1969. 

39. HITES, D.J. and MARCIL, P.7. , Determination of Upper 
Bounds for Problems in Plane Stress Using finite Ele- 
ment Techniques. Int. J. lech. Sci. 9, 5, 245-251, 

1967. 

40. ZIEMIEWIGZ, O.C. , and CHEUITG, I.K. , The Finite Ele- 
ment in Structiu?al and Continuum Mechanics, McGraw- 
Hill, Oxford, 1967. 

41. VEUBEEE, B.F. de and SiUCDER, G., ^In Equilibrium Model 
for Plate Bending Int. J, Solids and Struct. 4, 

447-468, 1968. 

42. MOKLEY, I.S.D. , The Triangular Equilibrium Element 
in the Solution of Plate Bending Problems. Aero.Qrt. 

19, 2, 149-169, 1968. 

43. HEREIiLillHT, L.R. , Finite Element Bending ilnalysis of 
Plates, J. Mech. Div. , ASCE, 93, EM 5, 1967. 

44. PR/iTO, C., A Mixed Finite Element Method for Thin 
Shell inalysis, Ph.D. Thesis, Dept, of C.E. , M.I.T., 

1968 . 

45. PAPENFUSS, S.W. , lateral Plate Deflection by Stiffness 
Matrix with Application to a Iferquee. M.S. Thesis, 

Dept, of C.E. , Univ. of Washington, Seattle, 1959. 

46. JOHHSOU, M.W. Jr. and MCLAY, R.W. , Convergence _ of the 
Finite Element Method in the Theory of Elasticity, 

JAM, E35, 2, 274-279, 1968. 

47. TONG, P, and PIAN, T.H.H., The Convergence of Finite 
Element Method in Solving linear Elastic Problems. 

Int. J. Solids and Struct. 3, 865-879, 1967. 

48. YI- YUAN YU, Generalized Hamilton’s Principle and 
Yariational Equation of Motion in Nonlinear Elasticity 
Theory, With Application to Plate Theory, J. Acou. 

Soc. Am., 36, 1, 111-120, 1964, 

49. DAYIDENKO, D.F. , An Approximate Solution of System of ^ 
Nonlinear Equations’, Ukrain. J. Math, 5, 196-206, 1953. 



91 


50. 


51. 


52. 


53. 


54. 


YiSOIiEV, M.I., Ihe Solution of Systems of Nonltoear 
Equations Toy a Method of Differentiation with Respect; 
to a Parani^Ster, U*S*S,R* Gomput* Mafh* and Ma-tli* Ptiys*^ 
PI, 4, 198-203, 1964. 


DAYIDEIKO, D.E. , An Application of the Method of 
Variation of Parameters to tTae Construction of Itera- 
tive Eormulas of Increased Accuracy for Numerical 
Solutions of Nonlinear Integral Equations, Soviet 
Mathematics, Doklady 6, 3, 702-706, 1965. 

SETIDR, A.V. and KAPOOR, M.P. , A Parametric Differentia- 
tion Method to Stru^ctural Optimization Problems, io be 
published in AIM Journal. 


JISCHEE, M.C. and BilRON, J.R. , Application of the 
Method of Parametric Differentiation to Radiative 
Gasdynamics. AIM 7, 7, 1326-1335, 1969. 


RJBBERI, P.E. and DANDAHL, M.I., Solutions of Transonic 
A.irflow Problem Through Parametric Differentiation. 

AIM, 5, 470-479, 1967. 



CHIPTER 3 


PLilE PROBIEM 


3.1 Introduction: 

Tro dimensional plane problem is one of the simpler 
cases encountered in continuum mechanics, though it involves, 
in general, high degree of nonlinearity for geometry as well 
as for physical material properties. This fact has made it 
rather difficult to obtain solution to any but the simplest 
problem by traditional methods. However, with the advent of 
high speed computers and development of finite element tech- 
nique, the problems can be attacked with comparative ease. 

G-eome trie ally and physically nonlinear problems par- 
ticularly with elastoplastic material have been studied by 
-1 

Pelippa, considering displacement models. Applications of 

this method to inelastic structures have been considered by 

2 3 4 *0 

Pope , Armen, Pifko and levine , Yamada et al. , Isakson" , 

6 

Marcal and King . Viscoelastic problems are investigated 

n o Q 10 

by Oden , ’Ahite , I'ried^, v/hile Oden and Zrass , Punjino 
and Ohsaka have studied thermal problems. Analysis of 
nonlinear geometrical problems are considered by Walker and 
Hall^^ and Oden'^'^, 

In this Chapter, the finite element equations for a 
mixed rectangular model in cartesian coordinate system will 
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be developed for plane problems having nonlinear geometrical 
and physical properties associated with coupled thermal effect. 
In section 2, some simplifying assumptions have been made while 
section 2 and 3 have been devoted respectively for the develop- 
ment of the finite element equations and nimierical solutions of 
some problems, 

3.2 Simplification: 

The basic equations which has been derived in Chapter 2 
3,re quite general. Various practical restraints inhibit their 
use for the solution of engineering problems. For example, 
due to lack of experimental data, the properties of a material 
implied for dipolar stresses are not knovm. Moreover, the 
nonlinea^rities are bound to maJ&e the formulation quite compli- 
cated which needs justification for practical utility. Most 
severe restriction comes from the availability of computer 
time and the expenses to be incurred for satisfying academic 
inquisitiveness. Considering the above points, it seems 
justified to simplify the problem having the following res- 
trictions. 

1 . The mntion has been assumed to be quasistatic and any 
term involving acceleration will be neglected. 

2. The dipolar stress field is nonexistent. 

3. It is assumed that the resulting deformations are 
accompanied by large rotations. This implies that the squares 
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of the rotations are comparable to the strains. Thus nonlinea- 
rities are due to the squares of rotations in the strain-dis- 
placement Illations. 

4. Generally, the effect of deformation on thermal state 
is small. Hence the nonlinearity of strains in the thermal 
equation will he neglected. 


5, The constitutive relation for heat flux vector shall he 
given hy generalized form of Fourier law, such that, 

= K^3(q) ^ 

in cartesian coordinate system. Here, comma denotes the partial 
differentiation with respect to the subscript following it. 


6. In the equation of theimal state, the variation of the 

constitutive coefficients (e.g,, ) 

will he neglected. 


Considering the above restrictions, the constitutive 
equations can he prescribed as follows, 

✓s 

i(Di) = (E^^, e, tP) 


A. 

0, tP) (3.2-1) 





. tP) ■ 

and 


= (q) 0^^ 


where , 

®ki 


5 
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In the above constitutive equations dependence of past his- 
tory has been tacitly assumed through the argument t^. 


3.3 Finite Element Derivation: 

For the construction of finite element models consider 
the plane continuum to be represented by R number of rectangu- 
lar elements in cartesian coordinate system, with a typical 
element as m, a typical node as p and the number of nodes, 

M, surrounding a node p as shown in Fig. 3-1. Fhe independent 
field variables are displacements U^, stresses = 1^3) 

and temperature 0, where i, 3 take values of 1, 2. Hence the 
number of degrees of freedom per node is six and dimension of 
the stiffness matrices per element will be 24 x 24. The dis- 
tribution of each of the field variables has been assumed linear 
in X and X direction. Hence, interpolating function for a 
typical field variable, say Y, having nodal values as Y^jY^, 

Y^ and Y^ can be expressed in the form. 


Y 


< 0 . 


0 , 


Y 


1 


0 , 





where, 


and 


0^ = (ab - bx"* - bX^ + x''x^)/ab 
Og = (bx'’ - X^X^)/ab 
O5 = X^'x^ / ab 

0^ = (aX^ - x''x^)/ab 


(5.3-1) 


(3.3-2) 
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1 2 

where X and X are local coordinates and a, b are typical 
dimensions of an element (Fig. 3-2). The nodal displacements, 
stresses and temperatures have been arranged in the sequence. 




= 

^12 

^21 

U22 

U31 U ^2 

^41 

CM 



= <T;j'' 

m 12 

1 

rn22 

1 

T^^ 

m12 

^2 • •• 

m12 

4 

rp 22 

4 

and 

{©yh ^ 

= < 


®3 

®4 




where 


and 


D = 

1,2 

; 6 = 1, 


, 4 


.T 


values of the corresponding quantities and i, j denotes coor- 
dinate directions in conventional way. Distributions of 
shall be taken same as those of T^^. To obtain, explicitly, 
the matrices in the governing equations, the following steps — 
are necessary. 


1. Substituting the distributions (3.3-1) in (2.6-2)^ and 

taking into accomit that the surface integral in (2.6-2)^ shall 
have no effect on the assembled structure because of continuity 
of stresses and displacements, can be obtained after inte- 

gration and has been given in Table (3-1). Similarly, the 

matrix C3„-, has been shovm in Table (3-2) . The vector 11 can 
PP a 

be easily evaluated for different cases of body force and sur- 
face load distributions. 

2, Evaluation of 12 q : 

P 

Simplifying the expression for 12^ in (2.6-6), the 

incremental 12 o can be written as, 

P 
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* f w m(3i)* dA 

= ( “up 

" 1 m 


(3.3-4) 


■ ^ BIT 


i(3i)* = 4- T^^>* 


+!!a-) 


.. - 8 U. . 

1,3 3»i 


1 -(a) ( !!s_ + ^ 

■" T * • ( W-" 


(3.7-5) 


It oa« te .horn, that for the nonltaear etrain e^reseioh 
( 3 . 2 - 1 )^ cliosen for this case, 


^ = 26’^, 6 \ 


(3.3-6) 


„,ere 6% is Kronec.er delta. Hence =» written ae, 

J(3i)* = cOi)(^^) Ej^ t 

/S 

,(3i)(ia) - sliiE , 5(3i)(ki) =4'" 

Where - 


and 


i(3i) = ^ 


( 31 ) 




(3.5-8) 


W 


Sow after substituting the expression for ) from 

( 5 . 3 . 7 ) in terms of nodal field variables in E^. (3.3-4) .Kp 

C3JI *b^ Wl*i'fe"fc^^ 3^ y 


E2;= ll'pal " ll-x' ^ 


^ ,1-1 I 1 I and I If, I are given in Tables 3-3 to 

where j f l» ' "^pa '"PY 


(3.3-9) 
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H /i 

3-5 respectively. Matrix I' will "be same as I ' only the 

poc pet 

coefficients c^3i)(kl) -^q t^q replaced by C^^i) 

convenience, linear and nonlinear parts of the strain expres- 

Jt n 

sions have been separated in j i and | respectively, 

3. Expressions for and 

The free energy expression has been chosen as, 


Po P 





e) 


therefore, 


% 

dJi^ 

dt 



(3.3-10) 

where 0^^ 

= 

dijj .d . 

p 0 M-r ’ 

ID 

= PQ 6A. . ^80^ 

1 J] 

(3.3-11) 

and 0® 

= 

^0 ^ 


(3.3-12) 


Since the effect of deformations on thermal state is 
generally assumed to be small, the nonlinear part of E. . may 

X J 

be neglected for thermal equation. Moreover, the motion has 
been assumed to be quasistatic, i.e., the effect of A. . is 

X J 

not prominent relative to others. Hence substituting in (2.6~7)^ 
the simplified expression, 




o dt ^ 8 e 


ID 


+ C? e 


(3.3-15) 


-^1 


where E. . is the linear part of E . . and performing necessary 
ID ID 

integration, 04 , can be written as, 

YY 


04 


YY 


Y J 

YY 


yr 


(5.3-14) 
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where, 


YY 


^^2 


Syimn, 


(3.3-15) 


K, 'K, 


S ^ 


T ,l = 

~YY 


Syram. 




(3.3-16) 


(3.3-17) 


. ^ - t 9 10 (3.3-18) 

Kp = ’ ’ ■■” 

XT xiT 1- 'hfl've heeri giver ir Tahles 5—7 sod 
ard matrices { J^} ard {K^} have ceen g 

3-8 respect ively . 

After some matrix manipulation, tbe Incremental form 


of C4 , can he given hy, 
YY 


OaI. S^,i, = 


(3.3-19) 


where , 


v'^ 

Y a’ 


Symm. 


T , 5) X T 

Jg J 3 

T X Y X 5 

^5 6 t (3.3-20) 

. T X 5? 
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m m 

rr Tr 


K 2 


Syram, 


T „ T { } ^5.3-21) 

^8 ^9 


4. Expression for L4 y : 

Since fhe constitutive equation for heat flux vector 


has been accepted as 


= z^He) . e,. 


(3.3-22) 


the expression for 14^ can be written by substituting, proper 


distribution and integrating (2. 6 - 7 ) 2 * Thus, 

= I K 


(3.3-23) 


where the matrix | R^^,thas been given in Table (3-9). Assuming 
now that is approximately constant for a small variation 

of 9, the incremental form of 14^ can be vn:itten as. 


^YY ^ ^ ®Y 'N ^ 


(3.3-24) 


5, Expression for 16 ‘ 


Noting, 




(3.3-25) 


the expression for 16 can be written from (3.6-7)^ in the 


form, 




(3.3-26) 
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where I S | has heen given in Table ( 5- 10). Again assuming that 

the contribution of will be negligible in the thermal 

eq.uation, 16* can be written as, 

Y 

' 5 . 3 - 27 ) 

Substitution of the above expressions in Eqs, (2.8-1) 
to (2 . 8-3) yields the following incremental forms of the finite 
element equations: 


M 


2-^ ( 1 
m=1 


{ T* } 

- '11? ) = 

0 

(3.3-28) 

M 

1 ' ! 
m=1 

“PP'I 

{ T* 

) -(I I?l + I 




- 

I I 

' V' - 'V 


j> > = 0 (3.3-29) 


and 


r { ! C4 yy. 

m=1 



These equations can be solved by the method outlined in 
Chapter 2. 
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21) -2a 2b - a b a 

_.2b -2a 2b - a ha 

2b - a 2b -2a b ^a 

-2b - a 2b -2a b 2a 

. b - a t -2a' 2b 2a 

_ b - a "b -2a 2b 2a 

. b -2a b - a 2b a 

_ b -2a b - a 2b a 


- b 2a 

- b 2a 

- b a 

- b a 

-2b a 
-2b a 

_2b 2 a 

-2b 2 a 


4 


Table 3-1: Matrix 01^^ . 


4 


2 1 

2 1 


1 2 

1 2 

1 

2 1 


4 


4 


1 


4 


2 4 

Table 5-2: Matrix G3^pi. 
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12 


_2'bc’^^-2ac^^ 
-2TDC^^-2ac^^ 
-21)0^ ^-2ac^^ 


-21dc 


11 

31 


ac 


-213 - ac 


-213 c 


21 


ac 


13 

33 

23 


-2ac"'^-r2l3C 
.2ac^^-2lDC 
-2ac^^-2l3C 
1 ? 

^2ac l3c 

32 130 


-2ac^ - 


-2ac 


22 


l3C 


l3C 


11 


ac 


13 


ac^^- 1)C 


13 

2bc^^- ac^^ 

2ac^^- 

33 

2bc^^- ac^^ 

2ae32_ 

’23 

2bc2^_ ^^23 

2ae22. 

>3 

2bc^^-2ac^^ 

2ac''^.^: 

33 

2bc^^-2ac^^ 

2ac52. 

33 

2bc^"’-2ac^^ 

. 22 
2ac - 

33 

bc^^-2ac^^ 

.12 

ac - 


- be 
_ be 


33 


be 


23 


.2bc 
.2b c 


13 

33 


_2bc 


,23 


2bc 


13 


. ^<,35 . 

_ bc2‘'_ ac2^ - 
_ bc’’’’-2ac''^ 


32 

32 


ae''*"- be 


33 

23 


bc^^-2ac^^ 
bc^^-2ac^^ 


ac - be 

ac^^-2bc''^ 

32 o-ho53 


bc^^-2ac^^ 

31 o„33 


32 


bc^'-2ac' 


ac-^^-2bc 
ac^^-2bc 


be 


11 


ac 


13 


33 

,23 

13 


ac^^-2bc 

32 


,31 .33 


"be "^3 ac 


_ ac'^'^-2bc^^ bc'^'- ac 


21 


33 


ac^^- be 
ac52. ^,33 

ac22_ 13^25 


bc^''+ ac^^ 
bc^^+ ac^^ 
bc2"*+ ac2^ 

bc^^+2ac^^ 

bc^^+2ac^^ 

bc2*'+2ac^^ 

2bc^^+2ac"’^ 

2bc^'’+2ac^^ 

2bc2‘*+2ac2^ 

2bc”’^+ ac"*^ 
2bc^''+ ac^^ 
2bc2‘^+ g_(,23 


bc^^+2ac*'^ 

33 
23 


ac'*^+ bc"^^ 
ac^^+ bc^^ 
ac22+ bc25 

ac'^2+2bc'‘5 
ac52+2be55 
ac224-2bc23 

2ac''^+2bc^^ 

2ac^^+2bc^^ 

2ac22+2bc2^ 

2ac'’2+ bc'*^ 

2ac52+ bc53 _2bc5^+2ac55 


. bc^^+2ac^^ 

- bc2‘^+2ac 

» bc”'^+ ac"*^ 
_ bc^'*+ ac^^ 
_ bc2'^+ ac^^ 

_2bc"’'’+ ac^^ 
-2bc^V ac^^ 
-2bc2''+ ^q2'5 

-2bc^''+2ac"*^ 


2ac22+ bc25 


.2bc2V2ac2^ 


_ ac”^^+2bc^^ 

. ae52-.2be5^ 

- ac^^+2bc2^ 

^ ac'*2+ bc"*^ 
_ ac^^+ 

_ ac22-. bc2^ 

^2ac^^+ bc^^ 
^2ac^^+ bc^^ 

-2ac22+ bc2^ 

-2ac^24.2bc"'^ 

-2ac^2+2bc^^ 

-2ac22+2bc2^ 


Hob at ions: 

0(11)(11) = c(^2)(l1) 

q(11)(12) _ q13. q(12)(12) 
^(11)(22) ^ 

TiBIE 3-3: Matrix 


^ ^31. (,(22) (11) = cSI 
^ (,33’. c(22)(12) ^ c^3 

^ (,32’. (,(22) (22) = 0^2 



d g 

d'g^ 

^1n.2 

d g 

A^cr^ 

-d g 


d^g’ 

a^g^ 

-d^g? 

a5g' 

d5g5 

a5g2 

-d^g’ 

4''g^ 

d^g"^ 

a^g® .. 

-a'g^ 

d2g5 

d^g^ 

d^g® 

-a^g'f 

a^g5 

d^g"^ 

a^gS 

-a’gT 

,19 

d g 

d g 

d g 

"I 

-d g 

a^gS 

d g 

^2^10 
d g 

^2^11 
-d g 

a5g9 

dV 

^3„10 
d g 

-a'g'1 

a'g'’ 

d^g''5 

d g 

-d^g'5 

a2g'5 

d2g'5 

a'g^^ 

-a^g'^ 

a5g^5 

a5g'5 

d^g'^ 

-a3gi5 


2 

-d g 

a’'g'^ 

A^J 
-d g 

105 

-d''g'*' 

-a^gS 

a^g^ 

^2„1 
-d g 

-a2g+ 

Id^g^ 

a9g+ 

a'^J 

-d'^g 

-a^g'^ 

^1 6 
-d g 

^1n.8 

d g 

-d g-^ 

-a'g® 

-d^g® 

d^gS 

-d^g® 

-d'g® 

-a^g® 

a’g® 

-a5g5 

-a^g® 

-digi° 

aigi2 

-d'gS 

-a^g'^ 

-a2gi° 

a^g"'® 

-a2g9 

-a2g« 

-a5gio 

aV" 

-a\9 

-d’g-'^ 

-d g 

d g 

-d'g''^ 

-d g 

-a^g'^ 

,2 16 
d g 

-a2g« 

^2J6 
-d g 

-d^g'^ 

^ 3^16 

d-^g 

-d3g15 

^3„16 

-d'^g 


Notations : 

^2 ^ C^12)(11) + c(-^2)(22) 

^3 = c(22)(11) + C^22)(22) 


“ ife ^cxN ^ ® 

{ g ^} has heen defined in Tahle 3—6 « 


TiBLE 3-4: 


Matrix 



where 
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Coii'td, . . 



105 



16) 
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J. - 


1 

-n 




_L 

72 


1 '5 '1 

-61)0 -6ac''^ 

/ 

1 

_3bc -2ac^ 

f 

1 51 

-be -ac-^ 


-6'bc^-6ac^ 


_3bc^_2ac^ 


, 3 ^2 
-be — ac 


6’bc’*-2ac^ 


3bc'’-2ac^ 


1 3 

be -ac 


; 6bc^-2ac’ , 

1 

X 2 

5 bc^- 2 ac 1 

• T - x; 

bc^-ac^ 


• ■ r 

abc^+aac^ 

■ Jg = 75 < 

bcV2ac^ 

J 3 " 7 ^ 

1 3 

be '+ac^ 


2'bc^+2ac^ 


bc^+2ac ^ 


3 2 

be^+ac^ 


_26c''+6ac^ 


1 3 

_ bc'+2ac^ 

/ 

-be^+ac^ 


■z 2 

^ _26c'^+6ac ^ 


_ bc^+2ac j 


r7 ^ 

^_bc^+ac ^ 

/ 

/ _2bc'^-3ac^ 


^ _6bc -2ac'^ 


1 ^ 
/_2bc - ac^ 

1 *2! 0 

-2bc^-3ac^ 


% 2 
-_6bc^-2ac^ 


-2bc^- ac 

1 *5 

1 3 

2bc - ac^ 


6bc"'-6ac^ 


21DC - ac-^ 

3 2 
2b c'^- ac 


3 2 

6bc^-6ac 

V. T - 

2bc -3ac I 

I 13 

2bc‘+ ac'^ 

^ J? = 72 ■ 

i 1^3 

2bc ‘+6ac-" 

[ ^6- 7? 

5 1 3 

2bc'+3ac^ 

3 2 

2bc"^+ ac 


3 2 

2 bc'^+ 6 ac 


2bc^+3ac 
i ^ 

_2bc’^+3ac^. 


_2bc”'+2ac^ 


_2bcV ac^' 

r\ . 

2 

_2bc^+5ac 

4> 

% 2 
_2bc^+2ac 

V 

/ 

^-2bc^+ 

ae^j 
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'■ -bo''- ao’ 

■5 2 

ac 

1 3 

be ac^ 




J7 72 




_ J_ 

72 




3 2 

bc^- ac 

1 3 

be + ac 

3 2 

bc''^+ ac 

1 3 

-be + ac 

3 2 

_bc^+ ac 


bc'^-2ac^ 

3 2 

_ bc^-2ac'" 

1 3 

be -2ac-^ 

7 2 

bc^-2ac^ 

3bc’'+2ac^ 

3 2 

3bc^+2ac^ 

_3bc^+2ac^ 

3 2 

_3bc^+2ac‘^ 
V 




Wobations : 


1 _ 0II 

c = Cq 




no 


72 \ 


r = 

10 72 


1 3 

_2bc -2ac^ 

3 2 

_2bc^-2ac 

2bc’^-6ac^ 

3 2 

2bc'^-6ac 

6bc"^+6ac^ 

6bc^+6ac^ 

-6bc"'+2ac^ 

_6bc^+2ac'^ 


^ .2bc^-6ac^ 

rz 2 

-2bc'^-6ac 

2bc''-2ac^ 

3 2 

2bc^-2ac^ 




6bcV2ac^ ^ 

3 2 

6bc^+2ac 

_6bc”*+6ac^ 

_6bc^+6ac^ 




= c 


22 


3 


.12 


Table 3-7: Matrices (p - s •••? 10) ^ 





f 



tA 

tr\ 

tA 





'44 


rO 

rQ 

rQ 

42 


C\J 



CM 


+ 

+ 

+ 

4* 


CM 

CM 

CM 

CM 


CH 

u 

4*i 

44 



03 

03 

03 




CM 

CM 


1 

1 

1 

1 


K\ 

tr\ 


tr\ 


%-r 

U 


44 


03 

03 

03 

o3 


CM 



CM 


+ 

+ 

4- 

4- 


T~ 

T-“ 

T- 

T~ 





44 




42 

42 




CM 

CM 


1 

r 

1 

1 


tr\ 


K\ 

tA 




<H 

44 


rQ 

rO 

rQ 

42 



CM 

CM 



-J- 

+ 

4* 

4 


CM 

CM 

CM 

CM 




4-i 

44 


o3 

o3 

o5 

03 




CM 

CM 

CM 





T— 

trs 

t<^ 

tc\ 

to 

' / 



44 

44 

Bi 

o3 

05 

05 

05 

rt3 


CM 

CM 



4- 


4 

4 

11 

T*“ 

T“ 

X”* 

XT" 


<H 


«H 

44 

tA 


rO 

42 

42 

44 



CM 

CM 


K\ 

trs 



« <K 



44 

!S^ 



PQ 

42 

44 

CM 


CM 

CM 

42 

CM 

t 

I 

1 

1 

' — ^ 

CM 

CM 

CM 

CM 

E4 



4^ 

44 

rd 

03 

03 

o3 

03 

CM 

CM 



{! 

^r^ 


IS^ 

K\ 

CM 



44 

44 

4^ 

03 

o3 

o3 

o3 



CM 

CM 



1 

1 

1 

! 

• rv 

T~ 

T— 

T“ 

T~ 




44 

44 

t— 


P 

42 

42 

T~ 

CM 

CM 



• E4 

EA 

tA 

tA 

tA 

rd 



44 

44 


rQ 


42 

42 


CM 



CM 

11 

1 

t 

1 

1 


CM 

CM 

CM 

CM 




44 

4^ 

T~ 

03 

o3 

o3 

o3 

44 

CM 

1 


I 

1 


tc\ 


£C\ 

tA 

«# 


.CH 

44 


m 

o3 

o3 

05 






CM 

o 

1 

! 

1 

1 


T” 

T" ■ 

■X — 

XT* 

■P 

M 


44 

44 

03 

rQ 

42 

42 

42 

P 

CM 

CM 



O 

! 

1 

I 




CO 


•H 

u 


I 


<D 

H 

. 


CJ 
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3.4 Numerical Examples and Comments: 

5.4.1: To show the applicability of the finite element model, 
three problems have been solved for which the typical arrange- 
ments of elements are shown in Eig. (3-2). Apart from these 
examples, various other simple problems whose results are 
known had been solved to check the reliability of the pro- 
gram. The errors were within acceptable limit. 

Example 1 ; The classical problem for an anisotropic half 
space with a concentrated load is studied and the stress 
distribution is shown in Eig. ( 3 - 3 ), The material for the 
half space is chosen as Topaz and the linear material coeffi- 

<44 

c lent s have been given by Hearmon The modified eoeffi- 
1 3 

cients for the two dimensional case are, 

q(11)(11) ^ ^^^24 ; q(11)(22) ^ g(22)(11)^ ^_23 

q(12)(12) ^ . q(22)(22) ^ ^^^2 

The accuracy of the solution has been compared through 
statical checks and the deviations have been found to be 
not more than 4 to 5 percent. Moreover, on the line of 
symmetry, the finite element solution compares satisfacu 
torily with that given by^ Lekhnitskii . 

Example 2; A viscoelastic cantilever beam has been analysed 
assuming the constitutive equations to be, 

( 5 , 11 ^ ^22) ^ (E^^, E 22 ) +0.5 (E^, Egg) 

t'' 2 = 0.5 E^g + 0.25 E^2 
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This problem has been treated both for linear as well 
as nonlinear' geometry and the *maximimi deflections at the free 
end and the deflection along the beam at time t = 1.5 have 
been shown in fig. (5-^). The stresses, obtained from the 
solution, have excellent agreement with the theoretical 
values. But, since they do not have any special signific^Oe 
for viscoelasticity (determinate problem) they have not been 
drawn here. 

Example 3: Second Danilovskaya problem has been selected to 
show the accuracy of the model for thermal analysis. The 
finite element solution and exact values as shown by Nickell 
and Sackman , have been compared in fig, (3-5)* for the 
comparison the following values are assumed, 

pQ = 0.008 , a = 15.3 x 10“^ 

P qC^= 0.96 K = 0.04 

^ = 42.0 X 10^ V‘ = 4.0 X 10^ 

where p^, a, 0^ , E, X and y are respectively mass density, 
linear thermal expansion coefficient, specific heat per unit 
mass, thermal conductivity and isothermal Dame's coefficients. 

The convective heat transfer across the surface is 
assumed to follow the relation, 

= h(T _ T^) 

Where direction i is normal to the boundary surface, h is the 
film constant which is assumed as 0.78, T is the temperature 
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of the surface and which is the ambient outside temperature, 
is taken as 300. The boundary surface temperature is pres- 
cribed to wary linearly from 300 'at t = 0 to 600 at t- = 

10 

1.45 X 10“ and then it is kept constant. 

Prom the fig. (3-5), it is seen that U 2 and T agree 
satisfactorily upto t = I .60 x 10 . After that they start 

deviating rather sharply. But even before t = 1.45x10“^® the 
accuracy for stress is rather poor, the maximum error being 
about 8,0 percent. Unfortunately, due to limitation of Com- 
puter time, no further refinement for time increment is possible 
to demonstrate the error characteristics. 

3 . 4.2 Comments: 

Prom a few problems solved in this chapter, it is not 
possible to give any general idea regarding the error propa- 
gation for the particuilar model chosen here. Hovrever, it has 
been found that better results can be achieved by decreasing 
the element size and step length at the expense of execution 
time. This has particuilarly been verified for the solution 
having short range duration. But for long range solutions, 
by decreasing the step lengths, the results come nearer to 
the actual values, but it Jias not been possible to take the 
step length so small that the propagation error is negligible 
even after an appreciable lapse of time. 
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Moreover, during tbe solutions of various small prob- 
lem^- for checking the program, it has teen found that fixed 
end boundary conditions yield less accurate results than a 
free or simple supported end. This is possibly because of 
linear interpolations adopted for the model. 

Continuum problems ane, as sucli^ quite complex and 
their detailed derivations are extremely complicated and 
massive. For this reason, the higher order interpolation 
laws become virtually prohibitive and it becomes mandatory 
to use linear law in spite of its disadvantages. However, 
the author feels that the results can be substantially 
improved, at least near the boundary gone, if some interpo- 
lation law can be used to impose such tha,t the adjacent 
nodes can be made to follow certain law so that the bounda.ry 
conditions can be satisfied more accurately. 







FIG. (3-4) STRESS DISTRIBUTION FOR HALF SPACE DUE TO 
A CONCENTRATED LOAD 
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CHAPITER 4 


AXISYlfflETRIG SOUPS 


4.1 Introduc-tion: 

Por deep sea exploration vessels and for aerospace 
structures, the analysis of complex axisymmetric bodies of 
arbitrary shape subjected to theimal and mechanical loads 
is of considerable interest. Because of the complex nature 
of the governing differential equations, the stress analyst 
must rely on numerical techniques. Although, the finite 
difference technique has been most popular, but it has its 
own limitations for nonhomogeneous materials and for arbitrary 
boundary conditions. Pinite element method has already been 
applied for the analysis of these type of structures under 

-j 

suitable simplifications. Rashid and Clough first applied 

* 

this technique considering axisymmetric torous element having 
triangular cross-section v/ith isotropic linear material. 

Wilson has extended it to anisotropic material subjected to 
arbitrary mechanical and thermal loadings using fourier ana- 
lysis, though for examples, he has mainly concentrated on 
axis 3 nJimetric loads. Recently Utku has given explicit ex- 
pressions of stiffness matrix for the triangular ring shaped 
element. It will be further extended here to accommodate 
nonlinear geometry and physical properties. 
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In the present chapter, the general governing equa- 
tions (1.10-1) to (1.10-3) will be simplified in cylindrical 
coordinate system (Pig. 4—1) for materials having constitu- 
tive equations as functions of deformation Velocity gradients 
and temperature. The temperature distribution has been assumed 
to be known as in the case of uncoupled thermoelastic prob- 
lems, though the equation for thermal state with mechanical 
coupling can be derived without any difficulty. The basic 
motivation is to reduce the system to a two dimensional case 
by uncoupling the circumferential direction (9-coordinate). 

This can be done by expanding the field variables (vshich are 
three displacements and known temperature distribution) in 
harmonic series in terms of 0, It will be seen that the most 
general anisotropic material which can be incorporated here 
shall have thirteen coefficients (see Section 108-109 of Ref- 
erence 4). Moreover, for nonlinear cases, only axisymmetric 
deformations can be considered. 

In Section two, the general equations are obtained 
for the nonlinear case. Section three is devoted for linear 
case, while in Section four nonlinear geometry with linear 
physical law is considered. Section five will provide a note 
for nonlinear material properties and in Section six, illus- 
trative examples will be given. 
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4.2 G-eneral Equations and Notations: 

Adopting single subscripted variable system for stres- 
ses and strains, and substituting from Eq. (1.10-2) in 
terms of in (1.10-1), tne three governing equations can 
be rewritten in physical components for cylindrical polar 


coordinates. Thus. 


{ ^^1,1 -^-F ^4,2 -^^6,3 ^ 


Ti _ T^ 


“ + ^ (?1 - f,))n,dv 


+ / (P^ - n^T^ - n^T^ _ n^Tg)U^ ds = 0 


s 

2T 

r (fj + _L []} + T + 

''^4,1 r ^'2,2 ^5,3 r 


+ Pq (5*2 - f.|))U2dv 


+ / (P^ - n^T^ - - n.^T^) ds = 0 


^ i h,2 h,3 * 4 ^ P" a' 


-b '"3 ■ •"1^'"3 


+ / (P^ - n^Tg - n 2 T^ - n^T^) ds = 0 

(4.2-1) 

where 1, 2 and 3 corresponds to directions r, 9 and z 
respectively/, P^, E^ and f^ are boundary body and inertial 
forces respectively, n^^ are direction cosines of the normal 

• I 

at a point on the boundary surface s, are corresponding 


velocities on the 

boundary 

s, the 

single 

subscripted 

are, 





A 

1 11 

m = 

-2 

T 

-22 

h 

= m 
-33 

A. 

np = m 

■^4 "'12 

"5 = 

K^ 
CM 
< E-i 

'^6 

A, 

= T.|5 


(4.2-2) 
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and ^ represents the partial derivative of T-j with respect 
to r, ill quantities in Eq. (4.2-1) are physical components. 
Excepting the surface integral portions, Eqs. (4.2-1) are 
well known equations of equilibrium and are given in Section 
59 of reference 4. 


The constitutive Eqs. ( 1 . 10-6) , can now be written as, 

.A A A, 


T 


T, 


T, 


T 


T 

5 


A 

where T. 
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T 

^1 

+ 

T 

1 



- V^)/t + I5 

° 1.5 


C\J 


A 

T 

^2 


2 + IJp/r +- ' 

A A 

2,5 


A 

T 
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A 

T 

3 

3,3 


^^6 “3,1 



A. 

h 


_L 

""rT*" 

A 

1 I, 


1,2 - VA 




4- 

A 

h' 


1 + 

)/r) + T5 tr^ 

/V 

,3 6 

^ 2,3 1 

"5 

-J' 

1 

2 

I 

+ T3 


,1 



-f 

A 

^5 

{ (tjg 


■ TT _} + T. 
3,3 6 



A 

T. 
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1 

2 

I T 
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A 

,3 ■" h ^ 3,2 

/r 




A 

T 

5 


2 - up/r H. 


pi 

9 * 

( 4 . 2 - 3 ) 

A 

T. 

1 


• • * 9 

^6’ 1 ’ * * 

• ^6» 

(4 

. 2-A) 


and G^, 6g are classical strain, Yg are 

Rivlin-Ericksen tensors in modified form and T is the tempe- 
rature. The classical nonlinear strains (1.10-9) are, 
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^ h . ( D 3 ^^) 2)/2 

62 = (1X2^2 ^ 

63 = 1X3^3 + ((XJ^^ 3 )^ + (^2,5)^ + (U3^3)^)/2 

26^= (U,,2- V/^-^^2,1 ^ (^1,1^^1,2-V 

+ U2^, (1X2^2 ^ ^1> + ^3,1 S,2)A 

^^5 ^ ^2,3 ^ ^3,2/^ '*' ^^1,3^^'''l,2 " ^2^ ^2,3^^2,2'^^1^ 

2^6 = ^1,5 ^ %,1 ^ ^1,1 ^1,5 ^ ''2,1 ^2,3 ^ ''5,1 ''5,3 

(4.2-5) 

and the Rivlin-Ericksen tensors ( 1.10-11) for cylindrical 
coordinate in physical form are, 

^ = ^1,1 ^2 = (°2,2 ^ =^3,3 

2^ = '“i,1-°2)A 2^5 = ^2,3-"%, 

(4.2-6) 

It may he noted that Y’s are the time derivatives of the 
linear parts of S*s, 


To derive the finite element equations for the axisymme- 
tric solid, let the whole region he subdivided into R number 
of axisymmetric torous (ring shaped) element. Let each element 
have the shape of a polygon in radial section ( fig. 4-2) and 
let h he the number of nodes (N > 3) per element. A typical 
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nodal point may be denoted by p and element by m, A parti- 
cular node p may be common to s number of surrounding ele- 
ment. Eor an element m, the nodes are designated as 1,2,...,!? 
and tbeir coordinates ( z^), (r^, z^) f %)• -Also 

let the displacements IJ^ (i = 1,2 and 3) be expressed in tems 
of tbeir nodal values at N number of nodes through an inter- 
polation la?;. Thus separating the variables, the displace- 
ments can be written in the form, 

tJ^(r, e, z, t) = \(l, Z) (8) 

DjCr, 0, z, t) = B^(t) njT, 2) (e) 

Uj(r , 9, z, t) = Djr, z) (e) (4.2-7) 

and T (r , e, z, t) = z) (0) 

where (0) = cos n9 (0) = sin n0 and the repeated 
indices denote summation. The interpolation function 
^m ~ have been asstimed same for all the 

field variables, \j^(t), C^^Ct) and 

the nodal va]..ues of three displacements and temperature of 
the node m, at time t and for the nth T'ourier expansion of 
the loading and the indeces are m=1,2,..,M'; n=0,1,..,'”. 

The distribution for T is not very important, because 
it has been assumed that T is a known quantity. However, since 
this formulation can be directly coupled with a derivation 
for theimal state, and the combined set may be treated for 
coupled thermal equations, the above mentioned distribution 
for temperature will be adopted. 
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4.3 Linear Case: 

Neglecting the nonlinear terms in (4.2-5), the linear 
strain expressions can be written as, 


1 


.0 


\in 1 ® n 

n ”n TT*^ 


-3 - 

mn m, 3 

n 

26° = 


n -h B D )/r«B 
urn mn 

26° = 

^®mn ^m,3 

-n C D /r) 

ion Hr n 

II 

OVD 

VD 

CM 

\,3 

+ c D J 

mn m,v n 


i2 

H' 


Por this simplified case, the constitutive relations (4.2-3) 
and (4.2-4) can now be expressed as, 

T? = C. . e? + g. . Y. + C.T i ,3 = 1,,.., 6 (4.3-2) 

where for linear material properties, the material coeffi- 
cients C,. and are constants and are given by 


and 


13 

8 

°i3 = rrj ; 


C. . = C. 
13 3 


8 T. 

rt r- ^ * 

-i3 8 Y ^ ’ 


C. . = C.. 
-13 "31 


C. 

1 


0 

3T 


From physical reasoning, it can be considered that, 

= 0, i = 4, 5 and 6 (4.3-4) 

Prom the strain expressions (4.3-1) and constitutive 
relations (4,3-2), it is evident that the stress components 
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will not be uncoupled imless, 

" ^2,1 = S3 " S3 = 0 D = 4, 5 (4.3-5) 



The constants C. will have similar restrictions. Hence, 
altogether, there are thirteen constants from static consi- 
deration, thirteen constants from velocity consideration and 
three constants from temperature consideration. 


Subjected to the above restrictions (4.5_5), the stress 
components in (4.3-2) can be subdivided as, 

(4.3-6) 




it + TT + T? 
1 1 i 


where, 


T = C . . e . ; TT = C . . Y . 

1 ID D 3 ''ID D 


(4.3-7) 


and T 


X _ 


G.T 


1 1 

How substituting strain expressions from (4.3-1) in (4.3-7)^j 
TT can be written as. 




i = 1,2,3 and 6 


and Tj = -iL " ^ 


i = 4, 5 


(4. 3-8) 


(4.3-9) 


where , 


x1 = C., D ^ + C-o D /r + C.c D a/2 
im i1 m, 1 i2 ts / 16 m,3' 


hm = hs’h “A 
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^im ®15 \,3 * ^16 

^lm=-°14“n= 'V2r 

= (°U - V") ^ °15 ®m,3)/2 

syzr 


(4.5-10) 


■^im 


- G.^ 

i5 


Similarly, can "be written as, 


= (?L V + ?L ?L °mnK 


and 


?: 


i = 1,2,5 and 6 

B + X^ G )H^ 
im mn im ron'^ n 

i = 4, 5 


(4.3-11) 


(4.3-12) 


where _X?^ can be obtained from X^^ in (4,3-10) by replacing 


im 


*^ik ~ik* example, 

x1 = 0.^ D ^ + C.„ D /r + C.. D ,/2 

-.im ...11 m,1 ...12 m^ i6 m,3^ 

Prom (4.5-7)^ and (4.2-7)^, P? can be written as, 

P? = X. E h"' 

1 im mn n 


where X. = C D 
im 1 m 


(4.5-13) 

(4.3-14) 

(4.3-15) 


i = 1,2 and 3 

Substituting the expressions for stress from (4.3-6), 
(4.3-8), (4.3-9), (4.3-11), (4.3-12) and (4.3-14) in (4.2-1), 
and using distributions (4,2-7), the first equilibrium equa- 
tions can be written in the approximate foim (noting n^ = 0, 
since it is a complete ring element), 
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O A ’ • 

+ ^?m,1 + * ^L,3 B, 


+ ^^1M,1 - 4m + Hm,3 * (^lV^mVr}0, 
+ ” V ^ ®!m “nm ^^1in,1 

^ t’'li.-^2mVrf V ^ '>o(^ - I ^ 


3 v3 


- 4 ^ 4 < ^q SB = 0 

-I 

where D^ corresponds to hoimdary distribution, A represent 
area of cross-section and 3,, the element boundary for 


1 2 

constant. Since H and H are orthogonal functions such 

n n ° 


that, 


2 tc -* M 2% o o 

/ 4 nUo^f h2 h 2 ae = ^ 6„p 


(4.3-17) 


the Eq. (4.3_15) will be completely uncoupled with respect 

to 0. The symbol Z" denotes the summation over all the 

elements foiming the entire axisymmetric solid. In similar 

manner the other two equations can be obtained and since equa- 

« * # 

tions are valid for arbitrary A , B and C , (which ai^ 

pq. pq pq 

not equal to zero), the determinate number of equations will 
be obtained and can be solved for a particular n. Dropping 
the subscripts n, the three equations can be written in the 
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concise matrix form, 


+ (I,}p > 


m pni 

0 


m 


v/here , 


{ 2 } 


i Vi 

] B ■ 
< mn 

C 


mn 


{ 


V. 


m 


K 


and 


mn 



~11 

^12 

~13 


{ k) = ^ 

pm 

521 

~22 

^23 



?31 

?32 

v„ 

. pm 


K. 


= ^ K 


K, 


pm 


U> = 1 


„ V Kr-/- are as follows: 

Tlie expressions for iv^^, ...» 


K. = TC . 

11 ■ 1 


1 

/ (n^ X^jjj + ) ®p 1 


( 4 . 3 - 18 ) 

( 4 . 3 - 19 ) 

( 4 . 3 - 20 ) 

( 4 . 3 - 21 ) 

( 4 . 3 - 22 ) 

( 4 . 3 - 23 ) 

i )/r)B dA 
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K ^2 = 5 ^ I {(Xlm ^ ^ 6 id ,3 


fX^ X' I'^r'lD' dA 


K^3 = It 1 


K . =-n% 1 J ( 
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^^3 em-^ p 
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■ K 33 = Tt 1 { ( - X 5 ^ 21 /r + X 3^^3 + Xg^ 
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3 m ,3 611 / P 
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Similarly, K . . can be obtained from K. . by replacing X by 
in (4.3-24). Tbe expressions concerning tbe temperature 
effect are, 


'1 = * I / DpdA - / d; ai| 


^2 =-n,| / X^yr dAl 


(4.3-25) 


^M,3 S 3 ^3a “p 41 

J-X ± 


Neglecting the inertia forces, the load vector is given by, 


h = U I>p 4A + ( 

a 4r ir 

S2 = 4 "2 “P 1 ^2 “P ' 

h = '{ So ^ “p 4^^ “p 4ll 


(4.3-26) 


where 



Fj) 

2it 

= f 

(P,, Pp) 

4 4e 

i = 1 and 3 

and 

F^) 

0 

271 

= f 

(Pg- Sg) 

Sp 40 

(4.3-27) 


o 


Now, by choosing convenient numerical quadratures over the 
area, expressions (4.3-24) to (4.3-26) can be evaluated and 
then applying suitable integration technique, the matrix 
differential equation (4.3_1S) can be solved. 


It may be noted here that all the loading cases cannot 
be solved easily by the chosen displacement foim (4.2-7). I'or 
example, the axisymmetric torsional (circumferential) load 
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needs special modifications. This type of problems can be solved 

easily if the distributions for and T are chosen by inter- 

changing and i.e., 

(r, e, 2,, t) = A^(t) ])^(r, z) H^Ce) 

Ug (r, e, 2, t) = b^(r, z) H^(9) , . 


and 

for t 

1 . 

2 . 


(r, e, 2 , t) = G^(t) b^Cr, z) h2(0) 
T (r,_B, z, t) = E„^^(t) D^(p,. z) h 2(0) 


(4.3-28) 




is distribution, the modifications are as follows: 

In Expressions (4.3_24) and (4.3_25) change n by -n. 
Replace Expressions (4,3_27) by, 

VI 27: O 

= I (I’i, Pi) H- de , i = 1 and 3 


and 


271 


( 4 . 3 - 29 ) 


(I®, pg) = / (Pj, Pg) ae 


The remaining portion of the derivation Vvall be unaltered. 


4.4 Honlinear G-eometiy: 

It can be seen that if the distribution (4.2-7) is 
substituted in the nonlinear strain expressions (4.2-5), the 
0-coordinate will not be uncoupled, unless it is assumed that 
nonlinearity in strains arises due to axisymmetric deforma- 
tions, or in other v/ords, the axisymmetric deformations are 
relatively much more prominent than those of nonaxisyimiietric 
nature. This tantamount to the restrictions that, 
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wliere q and n 7 ^ 0 and k = 0, 1, 2, 

and the left hand quantities can he neglected when compared 
v/ith the right hand quantities. 


Under this suposition, the nonlinear strains can 
he written in the form, 


£ _ gu 4 . g 

i i i 


i = 1, 2, 


) 6 


(4.4-2) 


where G? are linear strains (4.3-1) and G^ are nonlinear 
parts v/hich may he obtained hy substituting the distributions 
(4.2-7) in (4.2-5) with the restrictions (4.4_1). Thus the 
incremental expressions for the strains are, 

®* = ( ( V • 4o -^mo * (®pn- 

* S ®pn)) ^,1 


So = ^ ( -^ (A + n B ) + A (^C + n B* ) )I) D 

2 2 ^ ^o ^ pn pn mo ^ pn pn m p n 
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3 ^ pn ^mo 
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mo ■ pn 
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N i-.n •mr^ ■mr> -r^n pjj -^q jjjq -pjj / 


B - U ^ H 
ra,3 p,3 n 


•sr* 

2G 


4 P . ®pn^ ^m, 1 ^p ®pn ^p,1 ^m^ ■‘W 

- 'ip.o -'V ^ SpAVlV“H.'’p,4^ 


mo '"pn 


- n (C„. cL + G.„_ d. ) D„ ^ 


pn mo m, 1 


(4.4-3) 
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20* = l(-(nA +B)B_D+B B„) A* 

5 r ^ pn pn' iii,3 p pn p,3 ^ m< 

-A ( n A* B „ B + e ’ (B ,B _B „B )) 

no ^ pn m, :) p pn ^ m,3 p p,3 

-n (C C* + C G* ) B _ D i 

^ mo pn pn mo m,3 p n 
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"mo 


20 




(6 ■^* 


'‘"mo P 


* A + 6" G* C +0* C ) 
m mo mo pn pn mo 


(B ,B .+D „D J} 
p,3 m,1 m,3 p,1 

f ? 

where 6 =0 when n = 0 otherwise 6 = 1. 


Now neglecting any third power of nodal displacements, 
the incremental stresses can be written in the form 

Tf* = T°* 4- f* i = 1 ,..., 6 (4.4-4) 

where T°* is the linear part and are obtained from (4.3-6) 

^ 

and (4,3-7). fhe expressions for can be obtained from 
(4.2-3) by using the linear stresses (4.3_8), (4.3-9), (4.3-11), 
(4.3-12) and (4.3_14) and the nonlinear part of strains ( 4 , 4 - 3 ). 
Thus following the restrictions (4.4_l) and neglecting the 
terms which are more than second order in nodal displacements, 
the nonlinear part of the incremental stresses can be written 
in the form. 
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(4.4-5) 
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where the index 

: i varies from 1 to 6. The multiplier in the 

braces in (4*4- 

5) is when i = i, 2, 3 and 6 and 

'wfesn 

i =, 4 and 5-. 

The symbols used in (4-. 4-5) are as 

follows i 
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D Tq 

xl 4. 1 ..~ \ 1+1 ® 1 

^2p * ^2p r mo * ' r p ' 
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^6p 2 ^ ^^'l ^!ii,1 ’‘’ ^p ^m,3^ ^mo ^^p ^m,3 


Tn the expressions (x1..4_6) to (4.4-11)j 6' is zero 
whenever n = 0 otherwise it is one. The symbols and 
^im heen already defined in (4.5-10) and (4.3-15).T^i® 

expressions for 5*^ &^p ®^®-y obtained by considering 

only the curly bracket portion of the corresponding and 
G^ respectively and by replacing X^^ by ^ typical 
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example has heen given in(4.4-6). The expressions for 


IcTi Tr 

y. and I. are as given below, 
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The index i in (4.4-12) to (4.4-16) taloes the value 1,2,3 
and 6. Correspondingly for i = 4 and 5. The following ex- 
pressions are valid. 
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TTae stresses (4.4_5) together with T? from the lineax 
analysis can now he suhstituted in the Eqs. (4.2-1) and 
subsequently solved using quadratures for the area and line 
integrations using numerical integration technique for 
solving the resulting differential equations in terms of 
nodal field variables. The expressions, as such, are extre- 
mely massive. But a great deal of simplification will be 
achieved if the distributions (4,2-7) are assumed to be linear. 


It is worth mentioning here that this particular case 
is valid for the nonlinear stability analysis of an axisymme- 
tric solid subjected to an axisymmetric load. Snap through 
and bifurcation phenomena can also be studied using these ^ 
equations. 

4.5 Nonlinear Material Property: 

It can be clearly seen that non-axisymroetric deforma- 
tion cannot be uncoupled even for small deformations, when 
physical nonlinearity is accounted for. It is only possible 
for the case of axisymmetric deformations. Material nonlineari- 
ty will not pose any additional difficulty in this case, except 
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the coefficients C. 0. . and C. will no longer he constants, 

hut will he fimctions of C. Y - • and T subjected to the res- 

13 ^ D 

trictions ( 4 . 3 _ 5 ) at each point. 

4.6 Applications: 

4 . 6 . 1 : For the solution of prohlems which has been illustrated 
in this section, consider a triangular ring element with li- 
near distribution of displacements and without anj 

temperature effect, i.e., T(r, d, z, t) may he assumed to he 
constants. For this particular case U will he equal to 3 and 
for an element, the nodes will he 1 , 2 and 3 and their coor- 
dinates (r^, z^), ( ^ 2 * ^ 2 ^ ^ ^3’ ^3^ shown in Fig. 

(4-2). 

Now% since the displacements are assumed linea.r over 
a triangular element, the interpolation functions are given 

< ^2 > = <1 r z > I g j (4.6-1) 

21/3 2A/3 2A/3 

! S 1 = 2l ^2-^3 ^ 3-^1 ^1-^2 (4.6-2) 

r 3~^2 ^ 1-^3 ^ 2-^1 

and 2A = (rg-r^ ) (z^-z^ ) - (r^-r^) (zg-z^) (4,6-3) 

Using the interpolation functions ( 4 . 6 - 1 ) in Section 3 and 4 
different Cases can he solved through the use of computers. 
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Example 1 : 

Ixisjnnmetric deformation of a linear elastic cylinder 
has been analyzed and compared with the analytical results 
from the text of l.E.H. Love^. The end cross-sections of the 
cylinder has been prevented from longitudinal movements but 
are free in other directions. The radial displacement and 
stresses T^ and Tg agree excellently with the analytical results 
The longitudinal normal stress T^ and its mean value 1.66? has 
been found to be in good agreement. In all the cases the error 
is less than 1 percent. Eor plotting the stresses, the average 
value has been calculated between an element and its adjacent 
inverted one within each square unit. The results have been 
shown in Eig. (4-3). 

Example 2 ; 

The same cylinder but now in the foim of a thin disc 
has been analyzed considering the pressure variation as cos(n9) 
for different values of n (e.g. , n = 2, 4, 6 and 10). As a 
typical example the displacements and stress distributions 
fox n = 4 have been shown in Fig. (4-4) and (4-5). For ana- 
lytical results, a small program ha^ been made based on the 
derivation given in Timoshenko’s text^. For displacements 
the error never exceeds 0.5 percent, but for stresses the 
discrepencies may be as high as 3 percent particularly if the 
displacement has a very flat curve. This has been found 
for the stress T^. This has happened, possibily because of the 
low Values of the stress which, in turn, has caused ill condi- 
tioning. 
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Example 3: # 

The problem of a rotating anisotropic disc has been 

6 

solved and compared with analytical results given by Sam Tang 
in Eig. (4-6), The dimensions of the disc and the arrangement 
of elements are same as the previous example. Bfegnesium has 
been chosen for the anisotropic material and its elastic cons- 
tants are listed in Eig. (4-6). Here, it is only necessary 

2 

to put the body force equal to the inertia force r, where 
is the mass per unit volume of the disc and w qs the angular 
velocity. The results agree excellently. 

Example 4 : 

An anisotropic sphere made of Magnesium and subjected 
to internal and external shear load is shown in Eig. (4-7). 

This problem has been analyzed by Chen and is compared with 
the finite element solution. Since the material symmetry 
is in spherical polar coordinate, the elastic constants have 
to be transformed to cylindrical coordinate. This is done 
at the center of gravity of each element by a slight modifica- 
tion of the program. The shear loads at the inner and outer 
surfaces and the circumferential displacements at the angu- 
lar distance 67.5, 78.5 and 90.0 degrees from the apex have 
been plotted in the Fig. (4-7), 

All the above four problems ai^ linear and require only 
one step for solution. In multistep analysis, e.g., for non- 
linear or viscoelastic case, the available computer time 
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becomes tbe most stringent limitation. Por tbis reason only 

two problem has been solved for nonlinear and viscoelastic 
case. 

Example 5: 

To investigate the nonlinear aspect of the model, an 

annular plate has been analyzed. The dimensions of the plate, 

its material constants and arrangement of the elements are 

shown in Pigure (4-8). Initially the problem was solved by 

taking double the element sizes, shown in the figure. It has 

been observed that if the step for incremental load is reduced 

from 2 (which has been shown in Figure by the symbol X) to 1 , 

the solution deviates away from the firm line representing the 

8 

nonlinear solution for thin plate as given by Gordom . 

This is possibly because of compensation of errors due to dis- 
cretisation and large step integration which has resulted for 
the step increment 2. The variation of displacement with load 
for different incremental steps, clearly shows the importance 
of proper step sizes. However, with the step of 1, quite satis- 
factory results have been obtained upto the loading point 5, 
although it can not be stated very definitely because the 
analytical results are based on thin plate theory. 

Example 6: 

A viscoelastic cylinder stiffened by an external elastic 
encasing with ablating inner surface has been analyzed by the 
finite element method and compared with analytical solution 
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presented Idj B.C.Ting . ytie dimensions of the cylinder and 

the element suhdivit ions have been shown in Pig. (4-91). 

The Young's modulus and the Poisson’s i^tio of the outer 
elastic encasing are 3oe8pQ0-tiT5'e]_y -^o.xlO^ and 0,25. Now if 
a(t) is the inner radius at time t such that n(0)=5' ancl since 
the outer radius of innen sliell b is 10, the ablating rate 
from the reference 9 has been taken as, 

a^(t) = a^(0)/|-| j 

= 25/(l-3t/4) , t < 1 

and the pressure variatiojj with tespect to time has been 
assumed as, 

p(t) = (1 - 

The constitutive relation for the viscoelastic material has 
been prescribed as, 

^2^ ~ i ■*" €^idt' 

' ^ 0 12 

1 * * I 

+ / K(t - t') e,, e. j dt* 

0 id 

and ^ 

^5 = / d('fc-t')l€^+€„idt' 

"^0 i ^ 

where K(t) and J(t) for this problem are given by, 

K(t) = 82.0 9282.0 26 1 

and J(t) = 12472. 7 3094.06-'**'’^^’*^ 
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Due to ablation of the inner surface, at eveiy new 
integration step, the arrangement of the elements haire to 
be modified. This is done by the reducing the inner most 
column of elements and interpolating the previous values of 
elements and interpolating the previous values of strains 
at the center of gravity of new reduced elements. The inte- 
grations necessary for obtaining the material coefficients 
have been achieved by trapezoidal rule. It may be observed 
from the curve C of figure (4-9) that stresses and displace- 
ment at the interface are agreeable upto t '= 0.6. After 

this, the solution stairts deviating away from the actual 
solution as shown by firm lines in the figures. Also the 
stress distributions are quite satisfactory except near the 
interface. 

4.6.2: Comments: 

from the above analyses, it has been found that one 
step solution -is not at all a problem. But for multistep 
analysis and particularly for nonlinear geometry, the ele- 
ment configuration and the step sizes might be critical con- 
trolling factora. Since for the initial steps, a fairly con- 
servative step length produces quite satisfactory result, 
it is essential that for subsequent steps a thorough check 
has to be made. Obviously for a problem with known result 
it is possible, but for unknown problems checking the solu- 
tion by reducing the step sizes seems only possible. 
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Unfortunately it is a time consuming affair. However, 

10 

Davidenko Has sHown that the convergency can He mucH 
improved Hy employing more efficient integrating tecbniq^ue 
SUCH as Bunge-Kutta, HtHougH, it is a complicated situation. 
Hut in tHe long run it may He more economic. 
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FIG-(4-3) stresses AND DISPLACEMENT 
INFINITE CYLINDER (n--0) 
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CHAPTER 5 


TORSION INALYSIS 


5.1 Introduction: 

It has been recognised that stress analysis of a non- 
circular prismatic member of nonlinear material properties 
subjected to pure torsion presents foimidable mathematical 
difficiilties even in the simplest cases. Closed form solu- 
tions exist for shafts ha'ving only a certain class of stress 
strain law and mainly for circiilar and a few other cross 
sections. Energy approach coupled with mathematical or phy- 
sical discretization techniques may be applied for a more 
wider class of problems. In this chapter, the finite element 
technique has been applied for the analysis of torsion of 
prismatic members having arbitrary cross-section and varied 
nature of nonlinear material properties. 

In Section 2, a brief literature survey has been pre- 
sented. In Section 3, the complete form of three dimensional 
stress functions, introduced by Curtin, has been shown to be ■ 
identical to that of Prandtl for the torsion analysis. In the 
next Section, finite element approach has been employed for a 
prismatic member whose resultant stress and strain relation 
Can be prescribed by a continuous (e.g. arc sinh , Ramberg- 
Osgood etc.) functions. The cross-section of the bar has 
been divided into a finite number of triangular elements. 
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Within these elements, the stress fttnction 0, which satisfies 
the equilihriioii condition, is- expressed in terms of 0 alone 
01 ' intterms of 0 and its derivatives at the nodes, which in 
general, is termed as field variahles, The complementai^ 
energy is then expressed in terms of the unknown field varia- 
bles and using the principle of minimum complementary energy, 
a set of simultaneous nonlinear algebraic equations are 
obtained, lor a given value of the angle of twist 9, these 
nonlinear equations ai^ solved by step by step method in 
conjunction with an iterative scheme. In subsequent two 
Sections, examples are given for a triangular model with a 
linear distribution of 0 and then with a higher order dis- 

iq 

trihution, commonly known as KIT model in plate bending ana- 
lysis, with 0 Sind its two derivatives 0 and 0 at each nodes 

X y 

as field variables. It may be noted that the linear and higher 
order interpolation laws have their oito restrictions. The 
convergence of linear model is slow whereas the higher rrompa- 
tible model cannot be successfully employed for composite ma- 
terials. These short comings have lead to the development , in 
Section 7, of a mixed model having linear interpolation for both 
warping and stress function. The constitutive equations may 
he regarded as functionals with displacement and velocity gra- 
dients as independent variahles. The material symmetry will he 
restricted to a monoclinic system which has thirteen material 
coefficients, in general, for generalized Hooke ' s law which 
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results in three coefficients for plane stress torsion prob- 
lem. The last Section has been devoted for illustrative 
examples and discussions of results. 

5.2 Historical Review; 

Torsion of a cylindrical bar by end couple has a long 
history beginning with Columb and Cauchy . To Saint Yenant 
belongs the credit of bringing the problem of torsion and 
flexure of beams under a general theory by the inundation 
of the 'principle of elastic equivalence of statically equi- 
pollent system of loads.* Ror torsional analysis, he assumed 
that the state of strain consists of a simple twist combined 
with warping of the cross section so that the resultant stress 
at the ends are statically equivalent to a couple about the 
axis of the prism which ?/ill be equal to the applied torque. 
Later researches by Clebsh and Yoigt have resulted in consi- 
derable simplification of Saint Yenant* s analysis by the intro- 
duction of assumptions that there is no normal traction across 
any plane normal to the axis of prism and the state of stress 
is independent of the longitudinal coordinate along the prism. 

Basei on the above considerations, the small deflection theory 

1 

of elastic torsion can be derived easily . Adopting the car- 
tesian coordinates X , X and X^ such that X and X are in 
the plane of the cross section and X is parallel to the longi- 
tudinal axis of the prism and normal to the cross section 
(Pig. 5-1), the three displacements according to Saint Yenant ’s 
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and the ■boundary conditions, 

0 = 0g^ (Constants) (5.2-7) 

On the boundary Bi (i = 1, N) where Bi are simple non- 

intersecting boundaries of H fold connected region and G is 
the shear modulus of the material. From Eq.. (5.2-5) the 
applied torque T. can be shown to be, 

T = 2 / / 0 d x"' dX^ (5.2-8) 

the integration being extended over the cross-section. It is 
well knoTOi that stress function 0 is similar to the deflection 
of a membrane. This membrane analogy provide a convenient and 
intuitive visualization of the torsion stress function. 

With either of the two anproaches outlined above, the 
problem of torsion is reduced to standard problem in the theory 
of potentials in two dimensions. The most powerful tool for 
the analysis of potential problems comes from complex variable 
and associated approach which has been developed in the raonu- 
mental work by Muskhelishvili and a shorter account by Sokol- 
nikoff^. 

Eor the case of anisotropic material, it has been shown 

that a simple affine transformation reduces the problem similar 

to an isotropic case but with a different cross-section of the 
14 - 

prism ’ . 
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The prohlem of stress distrihution in a plastic cross 
section of a prismatic bar subjected to torsion is an inter- 
nally statically determinate problem. In the case of homo- 

* 5 

geneous isotropic rigid perfectly plastic material, ladai 
showed that the stress function 0 must satisfy, 

I grad 0 j - Is. (5.2-9) 

v;ithin a simply connected cross section, where k is the limit 
stress of the material together v/ith Eq. (5.2-7) on the bound- 
aries (where i = 1) and thus proposed his sand-hill analogy. 
Eor elastoplastic solution, combined membrane and sand-hill 
analogy may be employed. In 1941, Sadowskey^ extended the 
sand-hill analogy for multiply connected cross section. But 
as discussed by Prager and Hodge the combined membrane and 

sand-hill analogy is -valid ‘only if the plastic region remains 

8 

in the domain of influence of the boundary curves, Hodge has 
shown the plastic unloading (where stresses do not increase 
monotonically with angle of twist) aspect of multiply connected 
cross sections. This phenomena still requires further investi- 
gation to findout the detailed criteria for the unloading. 

Elasto-plastic stress analysis presents formidable 
mathematical difficulties even in very simple cases. Plastic 

Q 

zone may be evaluated by slip line^, but the elastic zone 
stress defies almost any kind of theoretical approach owing 
to the complex shape of the zone bouindary. For this reason. 
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practically all loiovm examples use inverse methods, where 
the boundary is derived from an assumed solution and an 

assumed elastoplastic separation line. Essentially the 

1 o 

sane indirect method has been discussed by Sokolovsky 
1 1 

ana Zener , Due to absence of any direct method, an appro- 
ximation formula for warping under large twist has been 
12 

given by Hodge , 

I>ue to these complexities, the ;Saint Tenant torsion 
problem acquires each year several new solutions by loiown 
methods. The conformal representation of assignable shapes 
has less significance, hecanse of the diversified nature of 
the practical problems and relative ease with which torsion 
problems can be solved, in the practical sense, by finite 
difference or finite element method. Because of its versa- 
tility and since irregular shapes, constitutive relations 
and material nonhomogeity do not pose any further difficulty, 
finite element method is more preferable. 

2.0 

Herimann first a.pplied finite element method for the 

torsion analysis of elastic prism with linear warping func- 

21 

tion over a triangular element. Shah et al. utilised the 
minimum complementary energy principle for the generation 
of finite element equation for nonlinear material with linear 
distribution of stress function over a triangular element. 
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5.3 Completeness of Prajidtl’s Stress Function: 

Consider a material body of N-fold connected region 
and its interior and boundaries denoted by D and Bi v/bere 
i = 1,2, I and each Bi represents a simple noninter- 

secting surface. In tbe absence of body forces and couple 
stresses, the Caucby*s laws of motion, i.e., tbe stress 
eg_uations of equilibrium for tbe continuous medium reduces 
to, 




mk km 


(5.3-1) 


where 60“^ is symmetric stress tensor and referred to 
rectangular cartesian coordinate system. A general solution 

of (5.3-1) is a representation of ?/itb sufficient functional 

2 

arbitrariness to span tbe class of all €0 in D that satis- 
fy (5-3.1). Ibe first stress function solution of (5.3-1) was 
given by Airy for two dimensional case and Maxwell (1870) and 
Morera (1892) for three dimensional case, Beltrami's gene- 
ralization”^^’ in tensorial form, may be read as, 

(5.5-2) 


m =6 0 A -A 

km kip mjq "^ijjpq ’ ij 




CoiTipleteness proofs for Beltrami's stress solution were given 

i R 

by many authors, e.g., Bom and Scbild ^ which all valid 
for simply connected region. Rieder has shown that re- 
presentation is incomplete and can utmost provide solutions 
of (5.2-1) which represent totall3’‘ self equilibrated force 
system on ea.ch surface Bi. Gurtin is the first to establish 
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a complete representation in terms of another symmetric second- 
order tensor fielf'^ and a 'vector field, v/hich admits the solu- 
tions of (5.3-1) and can be written ran, 

:2. 


ip 


-■] — 6-, . G A- - 

3cm kip EijCi i3,pq 


^3 = ^ 3 i = ’">1 = 0 


V‘='(IL „ + -H . .. 

k,m m,k' 3 , jkm 


(5.3-3) 


He has proved that this solution is complete in the sense 
that every (not necessarily totally self equilibrated) solu- 
tion of (5.3-1) may be represented in the form (5.5-3) for 
a If-fold connected region. The main a.ira of this section is 
to show that atleast for the Saint-Venent ' s classical torsion 
problem G-urtin's representation leads to Prandtl*s stress 
function. 


]From the classical Saint Venant assumption, the only 
nonvanishing shear stress components for a torsion problem 
are such that, 

T^^ = T^^ (x\x2) and (5.3-4) 

consequently the stress functions A., and H. in (5.5-3) will 
also be such that, 

" ^3 H 3 = H. (x\ X^) (5.3-5) 

Expanding (5.3-3) for T^^ and T^g in (5.3-4) with considera- 
tion of (5.3-5) will result in, 

= & 4. w fp ■ = + W 

^ 7 p ^ ^ f f '^’32 ^ H p 

and '^Ml "^22 ° (5.3-6) 
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with the hotaidary condition (on the assumption that the peri- 
pheral surface is stress free), 


where 


+ i! 

ds" dn 


‘23,1 " ^13,2 


'i' = H „ H 

3,11 3,22 


(5.3-7) 


(5.3-8) 


is the resultant fraction normal to the peripheral bound- 
ary surface and n and s are respectively the outward normal 
and tangential directions on the boundary surface at a point. 


IJow since is a plane harmonic function, which is 
apparent from Eq. (5.3-6)^, there exists a conjugate func- 
tion A which is such that 'i' + i A is a function of the com- 
1 ? 

plex variable X + iX . Hence from the Cauchy- Riemann rela- 
tion, if A is found, can be evaluated by means of the equa- 
tions, 


= A, 2 and - A,^ (5.3-9) 

Substituting (5.3-9) in (5.3-6), 

^ ^ ’ 2 *^32 ” (5.3-10) 

vfhile (5,3-6)^ will be atitomatically satisfied, The boundary 
condition (5.3-7) will yield, 


Tn = + A) = 0 (5.3-11) 

O.B 

In (5.3-10) and (5.3-11), if ^+A is replaced by ^ , then 
it shows that ^ is nothing but the stress function put 
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forward by Prandtl. Purt.ber, it is easy to sho'v that, 

*'*’22 " (5.3-12) 
Since Prandtl’s stress fimction bas to satisfy a Poisson’s 
equation as compatibility criteria, it is clear from (5.3-12) 
that A v/ill provide the homogeneous part of the solution, 
while ^ + A will supply the general solution. 

5.4 Finite Element Derivation Prom Minimum Complementary 
Energy: 

In this Section, .finite element equations will be 
derived through minimization of complementary energy in 
terms of nodal field variables for a triangular finite ele- 
ment model. It wil.l be assumed that the stress function 0 
?/ithin an element is expressible in terras of the nodal field 
variables through some interpolation law. The nodal field 
variables and che interpolation law will be left unspecified 
in this section and two special cases wall be derived in two 
subsequent sections. 

Similar to the cases of previous derivations, let the 
cross section of the member be discretised into R number of 
triangular elements having N number of nodal points. A typical 
nodal point ma.y be denoted by p and an element by m. A parti- 
cular node p will be common to M number of surrounding ele- 
ments. For an element m, the nodes are designated by 1,2 and 
3, their coordinates by (x!| , X^), (X^, X^) and (X^, x|) 
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(Pig. 5-2). Ill quantities which a.re asterisked represent 
nondiffientlonal form and their relations with dimensional 
quantities are given in Table (5_1), where ’a.’ is a typical 
dimension in the cross section. Stress fimction at any 
point within an element 'm* can be expressed in terms of 
the field vari.ables (nodal unknowns belonging to the element) 
by the relation, 

0* = <c > {h> (5.4-1) 

where, <?> is the interpolation function and ( b) is the 
chosen field variable vector appropriate for the element m. 
Differentiating the equation (5.4-1) with respect to X and 

p 

X , the stress matrix can be vritten down from Eq, (5.2-5). 

Thus, , 

[Ti^! 

{ a } = _L ' > = is!{b> (5.4-2) 

^o 1^25 1 

The constitutive equations relating the shear stress vs. shear 
strain is assumed to be monotonically increasing function and 

each volume element will be subjected to proportional loading 

1 8 

so that Henclcy's total strain energy "theory is applicable . 
The relation betv.'een resud-tant shear stress and shear strain 
can be expressed as, 

T = f( (5.4-3) 

where y^= (S^^)^ + ^®23^^ (^^ 5 )^+(T 23 ^^ 

and f is a single valued continuous function, and and 
may assume any value required for cu3n/e fitting and non- 
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diinensionalisation. From Hencky's total strain energy 


theory, 


^ 13/^13 ^ 23/^23 (5-4-4) 
From Eqs. ( 5 . 4 _ 3 ) ana ( 5 . 4 - 4 ), the strain expressions are 
given hy, 


0,2. f (5.4-5) 

and E 25 = 

where, 

The corresponding matrix relation in terms of nodal field 


variable (b) is 


{£ ) 


Pl3l_ "^0 

t fC 5 I T 


{b> 


(5.4-6) 


By definition, the complementary energy density for the 
element ’m’ becomes, 


A tt 


CO 


= T 


/ 

0 




d 


(5.4-7) 


Substituting Eqs. (5.4-2) and ( 5 . 4 - 6 ) into (5.4-7) and noting 
that the matrix I S 1 is invariant with respect to stresses 
and deformations, the complementary energy expression reduces 
to, 

AU = T Y ( ^^^ (h)^ ! G I d (h) ) (5.4-8) 

CO 0 0 Q 

where |C| = jSj ~ |S j (5.4-9) 

Using the resultant stress and strain relation, and Eq. (5.4-2), 
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it caB be easily shov.'n, 

dT* = {b}'^ let a {b>/ T* (5.4_10) 

Hence, the total complementary energy’’ for the m-th element 

can be obtained frora Eq. (5.4-8) -usine (5.4_10), Thus, 

I f ( '^*) j dx^ dy* (5.4-11) 
0 0 \ 0 ' 

where is the nondimensional area of the element and 1 

is the length of the privsmatic bar. The complementary work 
is given by , 

W = 219 :f f 0 dx* dy* = 219 <Z> ^b > (5.4-12) 


v/fcere < Z > 



<C > dx’‘ dy* 


( 5 . 4 - 13 ) 


Therefore, the complement. ary potential energy a for the 
overall member can be obtained from (5.4_11) and (5.4_12) 


after summing over all the R elements. 


■n: = 


R 


T Y_ la" 
o o 




( IS. 

r' 


'm 



-20* < Z > { b} ) 


Hence , 

f(T ) dv jdx dy 


Differentiating % with respect to the field variables for 
each node from 1 to H and equating each result to zero, R 
number of ma,trix equations will be obtained in N number of 
unknovm set of field v.ariobles. 


After differentiation, the finite element equation for 
the typical node p will be of the form. 
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s’ |Gj {b} = 29 * s'{Z} (5.4_14) 

where, 

}&l= ff. ICj dx* dy* (5.4-15)‘ 

K. "" 

and e’ is the s-ummation extending the M elements Vt^hich have 
the common node point p. 

It may be noted that for linearly elastic case f(T')/T 
assumes the constajit value v/hich is the inverse of shear 
modulus G and for this the resulting equations leads to ¥. 
number of linear algebraic matrix equations which can be 
solved directly. However, for the nonlinear form of f(T^ )/’^ 
standard step by step method can be adopted for the solution 
to be marched out from the unstressed state to a required 
value of the angle of twist. 


5.5 linear Model: 

Por linear model, each node has only one degree of 
freedom and the chosen field variable ^b^ is, 


01 


{b} = 


02 f 

j 

* 1 


0 


Ihe interpolation function < c > 


(5.5-1) 


can be derived to. 


<C > 


<^ > jg i 
* * 


* 


< 1 X -x^ y -y^ > 


(5.5-2) 


where < ? > 


(5.5-3) 
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_L 

2A 


and 


2A 


1 

0 

0 



K% 

1 

CM 

^’c 

* * 
y-j - 

* 

^1 - 

* i 

yg 1 

(5.5-4) 

1 41* •)f 

Ix^ — Xry 

* * 
x^ - x^ 

* 

Xg - 

* 

x^ 

1 

1 (7-^ - J^J 

- (x^ - 

X.,) (yg 

- y.,) 

(5.5-5) 


The mat 

rices 

I S j and 

jo jure respectively given by, 



|S1«: 

®31 

gjj*: 


(5.5-6) 



-§21 

-gjg -gj. 



and I 


O 





■gjl + 8,1 

®21®22 ®31®32 

®21^23 ®31%3 


10] .. 



S 22 ■" S 32 

§22^23 %2^33 

(5.5-7 


Symm. 


2 2 

§23 + §33 i 



The matrix < Z > can he evaluated from Bq. (5.4— 13) hy using 
(5.5-2), 

< Z > = <Z^ Zg Z^> 


(5.5-8) 


where Z^ = gg^ + + i 


I- §22 + ly §52 


X 


(5.5-9) 




^x ®23 V ®33 


and. 


ff^ (x" _ xp dx^ dy* 


// (y"" - yp dx^ dy* 

HI 
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Substituting these expression in Bq. (5.4-14) and differen- 
tia.ting with respect to 0^ v^here 0^^ is the value of the stress 
function at the typical node p, the required equations can be 
obtained. 

The principle advantage of linear model is that it is 
very simple and requires very limited memory space in computer. 
But the disadvantage is that, though it gives monotonic con- 
vergence, the rate of convergence is rather poor. However, 
some acceleration scheme, such Atkin's del-square or Richard- 
son's exinsipolation^^ may be applied to improve the conver- 
gency of the solution. 

5.6 Refined Model: 

In this model, the stress function 0 as well as its 
slopes 0^ and 0,, will also he continuous over the edges of an 
element. It is somewhat a-wlcvard to achieve slope compatibi- 
litj'' for a. triang’ular model; however this result has been 
achieved by dividing the element into three subelements. 

This procedure for deriving a slope compatible model was 
first developed by Clough and Tocher”''^ for plate bending 
analysis and is knora as HGl model. Steps leading to the 
relation corresponding to Eq. (5.5-2) are as follows. 

1) Each triangle can be subdivided into three subelements 

by choosing a point 'O' and Joining the vertices. The coordi- 
nate of 'O' ma.y be any where within the element but for con- 
venience it may be chosen as the center of gravity of the 

-Jf , 

triangle (x^, y^). 
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2. Transform tbe global coordinate (x^‘, y*) to (x*, y’) 

by shifting the origin to ’O' as shown in Fig. (5-3). 

5. The following steps have to be repeated for each 
subelement '1' (l = 1, 2 and 3) having nodes (m, n, 0) as in 
S'ig. (5-4). 


(a) Transform the coordinates from (x', y’) system to 
local system (x, y) such thct x and y axes are respectively 
parallel and noTOa.1 to the side mn. The transformation will 
be given by, 


X = 

x' cos 

8’ 

- y* 

sin 

9' 

II 

x' sin 

9' 

+ y’ 

cos 

0' 


( 5 . 6 - 1 ) 


where 

(b) 


where 

and 


tan e' = (y^ - yj / (x^^ - x^) 

Assume stress function 0* over the subelement 1 as, 

0* = <^ > r?- (5.6-2) 

< C > =<•] X y x^ 3^" y^ xy^ y^> (5.6-3) 

1 1 
<a > = < a2 • • - ®'q ^ 


Then the local field vaniahles may be given as, 


■ 0 


->r 


I 1 X y x^ xy y^ x^ xy^ 


* i 


( b) = 0,- |.= I 0 1 0 2x y 0 3x‘^ y^ 


T 

0 




^001 0 X 2y 0 2xy 3y 


{a^} 


(5.6-4) 


It is to be noted here that the stress function and its slopes 
will be uniquely defined along the edge nrn if the field 
variables { b) is prescribed only at the nodes m and n and 
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"bhis ensures -the slope compatibility. Now, calculate { b 1 at 
the five points m, n, 0, p and q. where p and q are the mid 
points of Or and en. 

c) Transform { b} to the global system {b) by the trans- 
formation, 

{ b} = 

d) At p and q normal slopes v/ill be required. Denoting 
p or q by r, the normal slope at r can be obtained as, 

^^y* at r ^ P» '1 

(5.6-6) 

where = tan""* ^-y ’ A' ) ] j, (5.6-7) 

and n^. is normal to the corresponding side at r. Hence, five 

vectors have been obtained for each subelement, e.g. , b^, b^, 
11 1 

b , b and b . The first three has three elements in each, 
while each of the last two has one element. Keeping this in 
mind, in general, the i-elation can be written, in the form, 

{ b^k } = 0^ {cc^} (iio suimriation over 1) (5*6-8) 

where 1 = 1,2 and 3; k = m, n, o, p and q. 

4) The vectors b^j^ are not all independent with respect 
to the overall element, because they have to satisfy inter- 
subelement compatibility which can be written in the following 


1 0 0 

0 cos 0’ sin 0’ 

0 -sin 0' cos 0’ 


{h} (5.6-5) 


manner, 
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that is { "b } represents the nodal vector of the field 
variables for the element, h^ is that for a particular node 
i (i = 1, 2 and 3) and 0 is a null vector. Performing some 
manipulation for reduction, it can he easily shown, 
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1 
a 

, 2 , 

3 
a 


H 


{■b } 


(5.6-13) 


where, 


H 


(5.6-14) 


(■“■11 “ -^^2 ^2 "^ 21 V 

jj - 1^2 -^^21 (■^ 11 "'' S 2 ^2 "^21 

Having obtained the relation (5.6-13), "tbe general pro- 
ceduie! c'utlined in Section 3 may be followed to establish the 
finite element equations. The matrix |sj may now be obtained 
from (5.6-2) such that, 

Is! = 1 S M H 1 (b } (5.6-15) 


where. 


-0 

0 


-1 0 -2x -y 0 3x^ 


0 


01 0 X 2y 0 2xy 5y" 


(5.6-16) 


linally the matrix will be given by, 


< C> = < C> I H ! (5.6-17) 

where <Z> is defined in (5.6-3). The main disadvantage of 
this model is the complica,ted expressions for integration to 
obtain ( G | and < Z >. The best way to eval'uate these matrices 
is to apply numerical integration technique over each subelement . 
lowever, if it is assumed that the elements are quite small in 
size and f(T*)/T* will be practically constant over an element, 
the computational economy can be immensely increased even at 
the expense of increase number of elements. 
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5.7 Continuum Approach: 

In this section, torsional formulation will he obtained 
for prismatic members of composite viscous solids having arbi- 
trary cross section. The material properties may vary across 
one or more surfaces as in reinforced concrete. In the previ- 
ous formulation v;ith linear stress function, this problem can 
be p-'lved with minor modifications, but in the neighbourhood 
region of an interface between any two materials, the conver- 
gence and reliability may be too poor. In the 2nd formulation, 
due to continuity of slopes across the edges of an element, 
the model will not represent an actual distribution even at 
the limiting case of infinitesimal element size, where the 

I 

element edge is also an interface of different materials, because' 
in the actual case there will be jumps in stresses across the 
boundary. Motivated due to these trammels, a third model has 
been developed where the warping function as well as stress 
function has been chosen linearly inside each element and the 
equations (1.10-1) and (1.10-2) in Chapter 1, have been applied ; 
for the derivation of finite element equations. | 

Assuming the angle of twist 9 is a given quantity in | 
cartesian coordinate system (x\x^,X^), the displacements are, 



- e x^ x^ 
e x"* x^ 

< 1 x"’ x^ > 


g 



(5.7-1) 
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and the stress function is given by, 


1 x"” x^ 


'1 

$1 


> I g I < 

^2 

> 



^3 




^ J 


i g ! 

7/ill be given by, 

" 24/3 

2A/3 


2V3 " 

y2 y2 

■^2 

,r2 ^2 


^1-4 

y 1 y1 
ir3 2 

x]-x^ 

4-^1 _ 


(5.7-2) 


_L 

2:1 


and 2A = (X^ _ x]) (x| - X^) - (X^ - x] ) (X^ - X^) 

and and are respectively the nodal 

values of and 0. 

The stresses, as derived from 0, are constants within 
an element, 


w 


= T 


13 


= - < g, 


31 ^32 "=33 


/ ^ 
*1 

' *2 f 


1 


'23 


< §21 §22 ^23 H 


5 

2 

$ 

V 5 


; 


(5.7-3) 

The strain components, derived from displacement distributions, 
are. 
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E 


13 


w. 


< §21 S 22 ^23 > 1 ^2 r 


-GX^ 


^ w- 
i 5 


Y 


'23 


^ §51 §52 %5 ^ ^ ^ 

w 


> + ex 


1 


X ^ 


(5.7-4) 


In r'vis. (5^7-3) and (5.7-4) g are elements of tbe matrix j g |. 

tj 

Since stresses from Eq. ( 5 * 7 - 3 ) are constants within an i 
element, the governing equation ( 1.10-1) reduces to the form, 

£1 / (P^-n^'T^) U^-f(P2-n^T2)U2+(P^-n^ ds | = Oj 

( 5 . 7 - 5 ) : 

I 

where, , P2 and P^ are applied forces over the houndary of the ■ 

12 1 

element, S in directions corresponding to X , X and X-^; n^, ng I 
and n^ aa:^ out ward normals; the integration is to he .performed , 
over the houndary of each element and the summation extends 
over all the elements comprising the asserahly. Since the stress i 
function satisfies equilihriun conditions it ’will he shown later ! 
that (5.7-5) is a redundant equation. 

■ ■ I 

Correspondingly, after minor modifications, Ici. ( 1,10-2) [ 

li 

can he written in the form, 

£{ / ( - 1|) Y ^ + ( "^ - T 2 ) X 2 i = 0 (5.7-6) 

A I 

where, x ^ and x ^ are stresses derived from the constitutive 
equations with assumed displacement distribution. The integration 
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extends over tlie cross sectional area of each element and sum- 

• 0 

mation extends over all the elements, /i^ain, since and 
are arbitrary strain rates and can be varied independently 
hence ^5.7-6) cai". be degenerated to, 

S I / - T^) dil i = 0 (5.7-7) 

ii. 

and 

s ! / (t^ - y dA j = 0 (5.7-8) 

^ 1 . *1. <C- 

ii. 

The constitutive equations for. the material in an ele- 
ment can be prescribe'! in the general form, 

.»s 

*^-1 ““ ” fc ^) 

Tp = ^ 2 ’ ^^2’ ”^1’ ”^2’ (5.7-9) 

As usual, they will be restricted by the imo dynamical restric- 
tion v;hich will not be discussed here. Since, these equations 
are nonlinear, parametric differentiation technique may be 


applied for the solution of problem. Applying this over the 
constitutive equations, (5.7-9) reduces to, 



‘ 1 

= Si ''l Sp. ^2 

* 

11 1 

+ 0^2 

* 

^2 


and 

'' * 
s 

^'21 ^22 ^2 

+ ^21 ^1 

^22 

* 

1^2 

(5.7-10) 

where , 

c. . 

ly 

9y" ’ *^i3 

tj 

ax . 

ji 

^3 



(5.7-11) 


C. . 
13 

- “jl “13 

= C , 
31 

i, 3 ' 

- 1, 

2 
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Prom Eqn. (5.7-4) 

* 



, the incremental strains are given hy, 

'< ^ S22 4 823 •'?> - 

'■* + 532 "2 S33 '"3) ■'■ (5.7-12) 


and time derivatives of Y.| and Yg axe similar to (5.7-12) 

•H* • • 

where and 6 are to he replaced by and 0. Now, since 
at each step Iqs. (5.7-5) (5.7-7) and (5.7-8) are to he 
satisfied, after application of parametric differentiation 
they will he laduced to the forms, 

E| / (P*-n5T*)U^+(P2-n3X2)U2+(P*_n^^*-n2'^2)U^ ds j = 0 

(5.7-13) 

Sj / (t*-T!i) Y^ dA 1 =0 (5.7-14) 


and 


2! / 




Yo dA 


= 0 


(5.7-15) 


Substituting in Bqs. (5.7-14) and (5.7-15), relations 
for r from (5.7-3), t ^ from (5.7-10) and expressions for Y^ 
siralnr to (5.7-12), the two equations yield, 

2 1 / {(§3^ + g32 ^2 S33 

+C1I (gg., w* + g22 w* + g23 w* - e*X^) 

* * * *,.1 \ 

^^31 ^1 ®32 %3 ^3 ® ^ ^ 

(§2^ + §22 ^2 + €23 W3 + 0 X'^) 

+^12 (g3., w* + g32 wj + §33 W 3 + 0 *x'^)>x 
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^21 ^22 ^2 ®23 S } <31 I =0 (5.7-16) 

^! / ^(§21 ^1 + 222 ^2 ^ S33 '^3) 

“^21^^21 + §22 ^2 ■'■ ^23 ^3 

-Cgglg^l w^ + Wg + +0 X ) 


• * * * p V 

-CgiCggi ^’1 + §22 ^2 ®23 "''^3 "® ^ ^ 


-€22^231 + §32 ^2 + §33 +e X^) } 

( S3 “I -j 232 ^'*^2 233 ^3 63 C } dl ] 


= 0 

(5.7-17) 


These equations have to he satisfied for a.rhitrary 
values cf and w^. Hence, after integration equating 

the coefficients of v/^ , Wg and w„ separately to zero will 
lend to six governing finite element equations corresponding 
to six unknowns per element. The equations can he written 

in n matrix form for the nodfe i where, i may he 1, 2 or 3; 

« 

~ ’fl ""’l 

S a’ ( |Ct|j^ v/g + !G-lj_ Wg + |H 

V.U 

2 3 

(5.7-18) 

w'hcre 1 is the area of the element. Each jG]^, Ig||^ and jHj^ 
are 2 x 3 matrices and the elements of the matrices are as 






1 

^2 5> 
*3 


) = 


2(1 . > 
1 


follows: 



18G 


hj 

^2i 


^11 %3 

"12) 

®23 

= g3j_ 

'S20 

^21 %3 

Cgg) 

H, . 

= iT^ . 



(5.7-19) 

'I3 

o2i 

^33 



H 

2d 

®3i 

^23 

where 

3=1,2 and 3. 


= §2^ A {(+ - C^2^) 9 + (C^ -,y - G 9 > 

Ig = S-^i a' U- 'Gg^y + ^22^^ ®* ^“*^21^’*’^22 ®* ^ 

A’y = / X^ rIA pjid a'x = / x'* dA (5.7-20) 

A A 

Tbe elements G.. can "be obtained from G. . by replacing 

Xj XJ 

0^2 rjid Ggg by G^g ^22 respectively. 

Considering 0 0, the remaining portion, i.e., the coef- 

ficients of e in (5.7-16) and (5.7-17) gives rise to a more 
concise form, 

s| / (t* - T*) X^ d.l = 0 



Consider, nov.^, the Eq. (5.7-13). Since the bar is pris- 
ma iic, the surface boimdary of each element, e, can be sub- 

, It 

divided into two parts, s the vertical surface and s the 
horizontal surfa.ce in the plane of cross-section, Hence, 


I 

for 0 , n^ = 0 

It 

and for s , n^ = Bg = 0 and = 


1 
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Concentrating on the surface s" and integrating over the 
whole cross-section, (5.7-13) will yield, 

I* = 2 S 1 / $* dA I (5.7-22) 

A 

¥/here T is the increment in toi«ional moment on the surface 
s*’. Again, isolating the surface s*, the equation reduces to, 

Sj P* - n^ T* - ngT* U5 ds 1 = 0 

(5.7-23) 

It may be noted here that P^ is not actually an applied force 
and has arisen because of discretization technique in finite 
element method. However, it can be evaluated from the displace- 
ment distributions through constitutive equations (5.7-10). 

If that is done, it can be easily shown by applying the Gauss 
theorem for transforming surface integral to volume integral 
that if equations (5.7-18) are satisfied, Eq. (5.7-23) will 
be automatically satisfied. In this system since, no known 
external forces are consideretj, the forces evaluated from 
displacements through constitutive equations may he considered 
as loading. With this consideration and by Eq. (5.7-13), ^d 
(5.7-22) it is not difficult to show that (5.7-21) are redundant 
equations and torsional moment can he calculated hy either of 
the relations, 

5* = 2 E j f dl j = Zj / (T*x'* - dA I 

A A 


(5.7-24) 
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The overall equations can be obtained by appropriately 
superposing Eq. (5.7-18) for the entire assemblage and this 
will give rise to simultaneous first order differential equa- 
tions in terras of incremental nodal values of ^ and w. This 
Can be -olved by any standard numerical integration scheme. 

5.8 numerical Results: 

Rig. (5-5) shows the arrangements of elements in the 

cross-sectiens considered for examples. The shear stress- 

shear strain relations considered for different cases are 

shown in Rig. (5-6). Therein, the experimental curves 1, 2 

and 3 are for circular cross-sections and 4 for rectangular 

23 

cross-section as reported by Buchanan . The curve 5 is 
given by the equation, 

y* = Y'3 sinhh (5.8-1) 

and the curve 6 is the classical elasto-plastic relation. 

The curves C^, and of Rig. (5-7), correspond- 
ing to the stress-strain relations 1, 2 and 3 of Rig. (5-6) 
shov/ the plots of torque vs. tv/ist in non-dimensional form 
for unit circular cross-sections. These are compared with 
the curves chov/n by Buchanan using finite difference 
method. In Rig. (5-8), the results for rectangular cross- 
section are obtained using curve 4 of Rig. (5-6). Herein, 
the finite element solutions are obtained using a linear model. 
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In Pig, (5-9) plots of T vs. 0 are shovm for the 
cases of imit circular and unit square cross-sections employ- 
ing linear model and compared with the results of Smith and 

1 8 

Sidehottom , for the stress stra.in relation given hy Eq, 
(5.8-1), It is to he noted that the results of Smith and 
Sidehottom are exact for the circular case, while for the 
rectangular cross-section energy method is employed. The same 
problems have been solved using refined model. The values 
obtained are practically same as given by Smith and Sidehottom. 
The maximum error for square section is about 0.2 percent where- 
as for circular section it is about 0,7 percent. 

The results for elasto-plastic analysis of circular 
and rectangular cross-sections using the present analysis 
are compared with the exact analysis, and are shown in Fig. 
( 5 - 10 ). 

Fig. (5-11) shows the plots of torque vs. twist for 
an 1-section corresponding to curves 5 and 6 of Fig. (5-6). 

This has been solved using the linear model only. 

Elasto-plastic torsion of a doubly connected square 
cylinder has been solved and is shown in Fig, (5-12). The 
dimensions of the cylinder are the same as those given by 
stout and Hodge so that the results can be compared. It is 
seen in Fig. (5-12a) that except at the sharp turning portion, 
representing the elasto-plastic transformation, the curve 
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p.grGor: excellently with the result given by Stout nn;' Hovlro, 
However, the naxirnum discrepancy betv/ecn values of 2 and 3 
for 9 is not pore than 3 percent. But, unfoidunately the 
deviation in stress distribution along the diagonal line 
as shown in Pig. (5-121)) is quite appreciable. This is because, 
to save conputer tine, this diagonal line had been taken as 
a boxmdary from synrietry considemtions. It has always been 
observed that error percentage on a boundary line is rather 
high coinpa.red to other inside points. Even so, the nature of 
the curves are fairly sane and unloading aspect fron elastic 
zone is clearly indicated at 8 equal to 7.06. 

It nay be noted that analysis of nultiply connected 
region poses no problem, Since at internal boundary contour, 
the values of 0 bas to be the sane everywhere. This is readily 
accomplished by a modification of the governing matrix equa^ 
tion. 
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CHiPTER 6 


PLilTE BENDING MililSIS 


6.1 Introduction: 

Structural plates of different shapes have a wide app- 
lication in the aerospace industries, Ihe finite element 
procedure is ideally suited for such general configurations 
since the plate may he approximated as a series of simple- 
shaped elements. Perhaps, the plate herding prohlem is the 
only field where the finite element method has heen most 
widely applied through various element models. Earlier 
attempts for the anal 3 rsis TOre confined to rectangular models 
and upto 1965, no satisfactory triangular model had heen evo- 
lved. failure of the previous investigations have heen attri- 
buted due to lack of interelement compatihility and absence 
of completeness of the assumed distribution field which has 
been discussed earlier. 

Notable contributors for rectangular displacement models 

1 2 *5 

are Adini and Clough , Bogner, fox and Schmit , Dea’k and Plan 

etc, Clough and felippa and Yeubeke-^^ have used quadrilateral- 
elements which have very good convergence property. Success- 
ful triangular displacement models have been evolved by Clough 
6 7 

and loci^r , Basel ey et, al, and others; whereas the equili- 
brium models have been illustrated by Yeubeke and Sanders 
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and Morley^, Assigned stress hybrid model has been developed 

1 0 1 1 

by Pian ^ and Severn and Taylor applied this to rectangular 

1 2 

as well as triangular models, Yamamoto investigated the 
behavior of assumed displacement hybrid models using rectan- 
gular elements. The question of convergency for a displace- 

1 '5 

ment model has been discussed by Waits et al. , whereas the 

necessary interelement continuity conditions and convergency 

14 

for various model have been given by Plan and Tong . 

For many engineering problems, a knowledge of the 
stress distribution is of as equal importance as the dis- 
placement field. Hence the finite element model should not 
only be capable of predicting displacements with reasonable 
accuracy, but also should he able to give a fairly accurate 
stress field. In some particular case, this nay be achieved 
by cboosing proper mesh sizes and adopting suitable averag- 
ing techniques. However, for many problems, this procedure 
may be impractical even with higher order displacement ex- 
pansions. 1 likely procedure to over come this situation is 
to adopt hybrid models^®’ Even in this approach, due to 
discontinuities at the interfaces, the fidelity of the solu- 
tion may he doubtful. Perhaps, mixed element model is best 
suited for this ptirpose, 

15 16 

Mixed element models have been studied by Herrmann * 

4 rj Pi 

and Yisser ' for plate bending problems and by Prato for 
general shells including the effects of shear deformation. 
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1 c 

Herrmann Has dewloped a triangular model with linear dis- 
tributions ;''-:r transverse displacement and tending moments. 

In another paper , he has taken linear displacements and 
constant tnoraents ¥?ithin each element and has continuity of 
normal moment at the interfaces. Even with this simple model 
he has Bhovm that the method yields good results. Visser'^'^ 
has adopted a parabolic vaiying displacement field with linear 
moment distribution over a triangular element using Herrmann’s 
variational formulation. It has to be mentioned here that the 
shear forces in these approaches have to be calculated numeri- 
cally and will be, at least, one order less acctirate than the 
moment values even if interpolation is employed. 

In this Chapter, a triangular mixed finite element 
model has been developed for plate bending problems in which 
shear effects are included. L:inear distribution for all 
variables (transverse deflection, rotations, bending and 
twisting moments and transverse shears) is assimed and the 
resulting matrix equation is obtained by the use of Reissner’s 
variational principle A detailed error analysis of the 
approximate equations is made and the convergence of the scheme 
is proved. This formulation has the advantage that (a) all 
quantities of interest are obtained directly; (b) error analy- 
sis is straight forward and hence corrections could be obtained; 
(c) extensici: of this formulation to shells of arbitrary shapes 
using curvilinear coordinates is feasible; (d) complications, 
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such as, gradual variations of thickness, temperature effects, 
initial stresses and dynamic formulation can he handled. The 
main disadvantage is that the total number of unknowns at a 
node point is higher than that of the existing formulations. 
This is partially compensated by the requirement of smaller 
number of triangles to achieve the same accuracy. 

It may be noted here that though the plate equations 
have been derived through Reissner’s energy approach, the 
governing equations which have been derived in Chapter 1 are 
equally applicable and when they will be simplified to this 
level, they will produce the same set of equations as those 
obtained through Reissner’s energy. Since, the Reissner’s 
approach is quite well known, it has been considered here. 

6.2 ¥ariational Formulation: 

For linear theory of moderately thick plates, the 
variational formulation considering transverse shear defor- 
mation is fairly well known^^. Hence without going into 
details, only the final form, relevant for subsequent deri- 
vations, would be presented. The notations used in the sequel 
have been illustrated in Fig. (6-1), The non-dimensionalised 
variational equation in condensed form is, 

6 Jr = // { 6 15^ (Bl5+ L)} dx dy _ / 6 F ds 
A ^ 

- / 6 ^ D ds = 


0 


( 6 . 2 - 1 ) 
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where B is an 8 x 8 matrix of differential operators given by 


. 2 - 2 ) 


the nondimensional state vector, 

# = < Pg w M^2 ^2 ^'2 (6.2-5) 

the nondimensional load vector, 

L® = < 1^2 - q 0 0 0 0 0 > (6.2-4) 

The relation between actual (primed) aiid nondimensional quan- 
tities are given by, 

X = x7h, y = y’A, = P^ §2 ^ ^2 

w *= w’/h 

= M*/Eh^, M2 = M^/Eh^, M^2 = (6.2-5) 

= V'^/Eh, Yg = ^ 2 ^^* 

= m^/Bh, nig = m|/Eh and q = q’/E 
where h is the thickness, E is Young's modulus and is 
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Poisson's ratio. It may be noted that the effect of normal 
stress is neglected here and hence the additional terms in 
L normally seen in the Reissner’s formulation are omitted, 

The surface integral in Eq, (6.2-1) extends o'^er the 
entire plate, the first line integral is over the boundary, 

s 

Q , where stresses are prescribed and the second line inte- 

s 

gral is over the boundary, d, where displacements are pres- 
cribed, The respective deformation and stress resultant 

s s 

vectors on o and d are given by, 

w > on ^ 

i)^ = < w - w* - p* > on h 

< “n - < “nt - ”^t 

Mt > on ^d 

where M_, 1 and are the nondimens ion al normal moment, 
n ’ n '^n 

shear and rotation along the edge respectively and the 
asterisked quantities are the corresponding prescribed non- 
dimensional values. The Euler Equations of (6.2-1) will 
furnish the first three equations as differential equations 
of equilibrium and the next five equations as force deforma- 
tion relations. The resulting eight equations can be written 
in matrix form as. 


B 1 + L 


0 


(6.2-7) 
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The independent vanishing of the line integrals will provide 
the required boundary conditions. 

”n = “S' “nt = “St \ = K “ 

and w = w*, = p* and p^ = p| on 

6.5 Finite Element Formulation: 

Within each element, linear distribution has been 
adopted for the approximate state variable vector, D. It 
may be noted that the total number of variables used here 
is larger than that required for a compatible displacement 
model. However, it has the advantage that, because of linear 
distribution, continuity of deflection, bending slopes and 
stresses will be obtained and due to complete pol 3 momial 

14 20 

distribution function, completeness and invariance criteria * 
will be automatically satisfied. However, there will be a 
limitation for this model, when the actual stress field is 
discontinuous, in indirect approach has been discussed to 
overcome this difficulty. 

To facilitate the choice of trial function, the plate 
is represented by R number of triangular elements having a 
total number of nodes equal to P. A typical node, denoted 
by p, will be common to M number of elements, as shown in 
Fig. ( 6 - 2 ) . For an element, the nodes are designated as 1,2 
and 3. Local coordinates and dimensions of the element have 
been shown in Fig. (6-3). A typical quantity, H, at any 


( 6 , 2 - 8 ) 
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1 ? : 

point (x , X ) inside an element, m, may be expressed in terms 
of its nodal values as, 


where 


t = 


2A 


21 = (a, 


H = 0^ t 


( 6 . 5 - 1 ) 

0 ^ = <1 xVa 

, x^/a > 

(6 

.3-2) 

h“ = < Hg 

H 3 > 

(6 

.3-5) 

(2/5)1 

(2/5)1 

(2/5)1 


( ^2 2 "” ^*5 2 ^ 

^ 32 ^ ^ 

— 3 . 22 / S’ 

(6.3 

; ( - 1 — 3-2 j )/s 


ag^/ a 


21 ^32 ^31 ^ 22 ^^^ ’ 




and ’a’ is a typical dimension of the plate. Following 
Eq. (6.5-1) the distribution of the state vector in terms of 
its nodal value may be given as, 

p = fl AD® (6.5-5) 

Where iF = < ^ Y^> (6.5-6) 

$ = Pia 0 ^ 0 ^ 0^ 0^ 0^ 0^ 0^J 

(6.5-7) 

T = Pia f*t t t t t t t t } 

(6.5-8) 

P®^ = < P^ p| d| > (6.5-9) 

where is D at the node, i. Matrix - is a 24 3C 24. boolean 
transformtion matrix to rearrange the elements of the vector. 

< Pl1 P-12 Pi3 P 21 ^22 P 25 **• ^15 ^21 ''^22 "^25 > ^ 
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In the above, the first subscript denotes coordinate 
direction and the second one corresponds to the node number 
of the element. It should be noted that D and B are not 
identical. Vector B represents the exact solution surfaces 
and is the solution of the differential equation (6,2-7) with 
boundary conditions (6.2-8); whereas, B is the approximate 
solution of the finite element method and is valid only within 
the m-th element. Hence the exact state vector ^ for element 
m, may be written as, 

B= $TfiBp + e“ (6.3-10) 

where, 6®“ is the associated error vector for the element m. 
Substituting B from Eq. (6.3-10) in the surface integral 
portion of Eq. (6.2-1), we get for the element m, the expres- 
sion, 

T ■ T — 

6 B® (|K® B^ + j + iP) + // j BB+ I i dx^ dx^ 

where = qV (ff B 

= flV iff L 

E® = Q ^ B e“ ax'® dx^ ) (6.3-14) 

and is the area of the m^th element. The 24 x 24 matrix, 
K®, and, for linear distribution of loaa, the exact expres- 
sion for may be calculated niauBrically. The expression 


dx^ dx^ ) T 


(6.3-11) 

( 6 . 3 - 12 ) 


dx^ dx^ ) 


(6.3-13) 
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(6.5-11) can be written in a slightly different form as. 


£ 

< 

> 

1 s 


E??. 

D? + + 1 

»2,3 


1 

3=1 

.2,3 

1-3 

3 1 


-1- 

// 66“^ 

I 115 + 

L 1 dx”* dx' 
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:gm 

= 

< E? 

1 

e“ 


> 



(6.5-15) 


( 6 . 5 - 16 ) 

TO. ni 

and partitioning of K and L matrices have been given in 
Tables (6-1) to (6-5). 


Performing snimnation over all the elements after 
appropriate superposition and noting that there are only M 
elements surrounding a typical p-th node, the expression 
(6.3-15) gives, 



T 

+ Z // 6€“ { 1 5 + L j dx^ dx^ 

m= 1 , . . R 

(6.3-17) 

where the i-th node of the element m is p. 


The independent vanishing of the coefficient of 6D 


T 


in (6,3-17) gives, 


Z 

m=1 , , 


.1 


C £ E? . D? + L® + E?) = 0 (6.3-18) 

3=1. ,2, 3 ^3 3 i i 


and Tanisblng of 66^ gi'^es Eq. (6,2-7)* 
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Consider the line integral in Eq, (6,2-1). Designate 

a typical quantity by H the conjugate quantity (with res- 

pect to work) by Y and associated error terms by end 6'^ . 

Integrating from node p to p + i along the boundary whose 

•1 

outer normal subtends an angle 6 with x axis, as shovm in 
Dig. (6-1), a typical line integra.1 may be written as, 

/ 1 H* - H - i 6(Y+ e'^)ds (6.3-19) 

s 

where H* is the specified value of H along s. Since H and 
Y are assumed to be linear functions, the expression 
(6,3-19) may be written as, 

I / (H* - H - s _ gH, _ ^) ds I 

s P ®P =p p 

+ I / (H* - H B _ eH) ^ as 1 6^ 

^ p p 

+ / 5 (H* - H) ds (6.3-20) 

s 

where s is the length of the path of integration from p to 

P 

p + 1 and the barred quantity is the corresponding exact val'-.e 
on bounda,ry s. Since the distribution of H* is known, inte- 
grations for the coefficient of 6 y ^ and may be easily 

performed. It is to be noted that Independent vanishing of 
the coefficient for will give the exact bomdary condi- 

tion (6,2-3), Using standard transformation relations to 
convert the normal and tangential components (n - t coordi- 
nate xsystem) to usual quantities, in (x , x system) the 
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coefficient of 6Y^ may be appropriately added to Eq. (6,3-18) 
to impose necessary boimdary constraints. 

Neglecting the error teims, the reduced matrix equation 
for the entire assemblage of elements which is of a banded 
nature, may be written as, 

K B + B =0 (6.3-21) 


where, K is the overall constrained ’stiffness’ matrix, 


ffi rp m 

B = < iq Eg 


... >, 
p » 


.T _ . jT 

j - < Ig 


T 

> 


( 6 . 3 - 22 ) 


Solving Eq. (6.3-21) by standard algorithms, such as the 
Gauss elimination procedure, the approximate vector 1) may 
be obtained. 


6 , 4 Erro r inalys is : 

Accuracy and convergence study of a proposed finite 
element model is generally based on numerical evaluation 
of some problems and comparing the solution with the available 
analytical results. Such a procedure may be valuable in pro- 
viding a qualitative character of different models, but is 
basically deficient in general n^thematical justification, 

J.l. Wals et. al, have investigated convergence properties 
of several finite element displacement models based on classi- 
cal order of error analysis. 
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In the following, convergence and error propeirties of 
the proposed model are studied to obtain the p.pproxiraate 
discretization error involved in the analysis. The discreti- 
zation error is due to the omission of the term iff in Eq, (6.3~1S), 
which in turn, is a consequence of the approximation for the 
distribution of the state variables. 


By Taylor series expansion, the functional value of 
D at any point a (a = 1, 2, 3) can be expressed in terms of 
5 and its derivatives at point i(i =1,2, and 3). Subtracting 
appropriate error terms at the nodes, the approximate state 
vector D can be written as, 


D 


a 


= E. + e. - e 

1 1 a 


+ SZ 
3 ,k=1 ,2 


3 = 1,2 

.f ( 


ai ( N 

^3 ^ 3x3 




axS ex’" 


)• + 


(6.4-1) 


where x 


ai 


3 


a • - a. . and 6 . is discretization error at the 
a3 13 1 


node i. Concentrating the attention on node p (recalling p 
is the i-th node of element m), substituting (6,4-1) in 
Eq, (6.3-18) and since the equation, 

(K®. Dff + 1®) = 0 


S 

m=1 , , 


■ Z 

3=1,2, 5 


ID D 


has been already solved for approximte solution, it can be 
shown that. 


S Z 

m=1 , . .M 3 = 1 , 2, 3 
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-m 


i3 


(e, - e.) + Eff + 


a=1,2,3 3=1,2 


ae 


-r "la ' >1 * 


Z A. 

3=2,3..“ ^ 


0 ( 6 . 4 - 2 ) 
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where A. ’s are given by, 

tl 


Ag 

. 1 

z 

(xi 

ES X. X- 

£. • 

a= 1 , 2 , 3 

3,k=1,2 3 

b = 

. 1 
3! 

E 

EEE x?^ 


cx=1,2,3 

3, k, 1=1, 2 5 


^ dx3 dx^ i 

(6.4-3) 

ai- ai / 9^ v 

k 1 ia ax^dx^9x^ ^ 


and other A’s can be similarly written. For approximate error 
analysis, the vector may be approximated in a linear form 
similar to that of the state variable vector D. Following this 
approximation, Eqs. (6,3-14) and (6.3-16) will give, 

E? = e® + (6.4-4) 

Substituting (6.4-4) in Eq. (6.4-2) and differentiating the 
term (dS/dx^)^ for linear distribution, it may be shown that, 


m = 1 , . . . M 


e® + e“ + K®, e? + 

i1 1 xd d 13 3 


4. j = 0 


3=2, ..a 3 

(6.4-5) 


Eq_, (6.4-5) can be designated aS the unconstrained equation 
for the error vector at the node p, and the expression in 
the braces as pseudo load. These can be calculated approxi- 
mately from the known value of b by any numerical technique, 
Consider now the part of the expression (6«'5-19) involving 
6 y p. Following the procedum adopted for the siarface intfe-, 
gral, the expression containing 6 1 in Eq. (6.3-19) can be 


written as, 
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®p ^ ^ ~~^ ^ ®p ^ * * ■ S (6.4-6) 

where I is the exact expression for H*-H on s. Igain, for 

approximate solution of error, talcing linear distribution 
we get. 


I + (i)p + ^ ^ --■ I (6.4-7) 

Since the derivatives of 5 can he calculated approximately 
from the previous solution, this expression (after applying 
proper transformation to (x% x^) coordinate system) has to 
he added for appropriate boundary to the Eq. (6,4-5) as cons- 
traint. Hence the resulting reduced error equation can he 
solved. 


To show the characteristics of the model more clearly 
all a^-s's Eq. (6.4-5) may he divided by a scale factor ^ 
where, 


a . • — A a • ■ 

13 3.3 

Such that when A tends to zero, all the M number of ^odes 
surrounding a t 3 rpical node p, tend to shrink to p, keeping 
the relative dimensions of each triangular ‘‘«le®^f '(a. .) to 


he the same. It can he easily shown that when ^ - 0, Eq. 
(6.4-5) is of the form, 


f 
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Lim ff I ('B e) + X{ } + A^'{ } + ... I dil = 0 

X - 0 A ’ 

( 6 . 4 - 8 ) 

It is seen frori (6.4-3) that expressions in the braces always 
remain bounded as ^ tends to zero and hence Eq. (6.4-8) will 
yield trivial solution 8=0. (Actually, it v/ill be only 
Solution from the uniqueness consideration of elasticity). 
Convergence to the true solution is thus assured. Monotonicity 
will be achieved if the lower order derivatives in the pseudo 
loading terms retain same sign in the region boimded by the M 
elements as the element size ( ^ ) is reduced. 

It can be shown that for regular element configuration, 
the coefficient imiltiplied by ^ in the load term is identically 
zero. Thus, the solution converges quadrat ically when elements 
are regular. Also, if the state variable at the nodes of an 
element are expressed in terms of variables at the centroid 
of the element by Taylor's series expansion, then, it is seen 
that the coefficient of ^ again vanishes, thereby indicating 
that the solution at the centroid of the element is better 
behaved and converges quadrat ically. This result has been, 
demonstrated by several authors. 

Finally, near the boundary region, the solution con- 
verges linearly because of the presence of the linear teim 
in ^ , 

Similar conclusions can be drawn for the expression 
(6.4-7) which sho?.'s the error coninitted in discretizing the 
state vector at the boundary of the plate. 
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6.5 Discontinuity of Stress Field: 

Discontinuity of stress field in the elastic ra,nge nay 
result from the following situations, 

(a) Where there is an abrupt change in the thickness of 
the plate. 

(b) When the plate is made up of different naterials. 

Across the discontinuities, the normal moment and shear 
are continuous, whereas the tangential components of the stress- 
es are discontinuous. Independent tangential variables nay be 
assumed on either side of the discontinuity and the correspond- 
ing equations may be obtained by the independent vanishing of 
the expressions multiplying the variations of these tangential 
values. The result is that the number of unknowns and the 
corresponding number of equations increase at such disconti- 
nuities, llso, for such cpses, the band width of the result- 
ing simultaneous equations Eq, (6.5-21) would increase. 

6.6 N-umerical Results and Conclusions: 

6.6.1 lumerical Results: 

To show numerical convergence and error properties 
several problems have been solved. 

Fig. (6-5) shows deflectior w', moment 
M’ , and shear, for a simply supported square plate having 
2, 3» 4 and 6 divisions in a quarter of a plate as shovm in 
Fig, (6-4) (Example 1). It is seen that, inspite of the large 
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clui.icnt 5l.~,es, the results are \’yell behaved and close to the 
correct s:lution. Monotonicity of results is clearly seen. 

Fig. (6-6) shows deflection, w’, and moment, M’, for 

a rectangiilar plate with opposite edges simply supported, 

one edge fixed and the other edge free (Example 2), Four 

sub-divisions have been taken for one-hp.lf of the plate. This 

16 16 

example has been taken from Herrmann ’ . The solutions are 

compared with the series solution given in the references 
15 and 16. Except near the fixed edges, the finite element 
solution shows excellent agreement with the series solution. 

Fig, (6-7) shows deflection w*, moments , 1^, and 

shear for a square plate with two opposite edges clamiped 

and the other two simply supported for thickness to span 

ratio 0.01 and 0 . 1 . Except the moHKjnts near the fixed edge, ^ 

other quantities are very well behaved and quite close to 

the actual solution. Exact solutions for deflections, moments 

21 

and shears have been drawn for the thin plate . In the finite 
element solutions, differences of moments and shears between 
the thick and thin plates have been obtained between 2 to 3 
percent except, only, in the case of shear Tg slong y = a/2, 
where the maximum deviation is about 11 percent. Exact value 
of the maximum deflection and moment for a thick simply sup- 
ported plate has been taken from the reference 22. 

Lastly, error corrections have been applied using 
Eq, (6.4-3) io the finite element solutions and the results 
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have heen shovvn in Table (6-5). Since the error at any node 
has been asstcicd to be the maximum of the corresponding error 
anywhere in the plate, the wide deviation of the corrected 
Value seems tc be Justified. Nevertheless, the corrections 
furnish a bound which gives an approximate idea for the loca- 
tion of the exact solution. 

6.6,2 Conclusions: 

i. mixed finite element model for plate problems has 
been developed. This model takes into account the effects due 
to shear deformation. Because of the linear distribution of 
the state variables, inter element compatibility is completely 
satisfied. Equilibrium and compatibility conditions within 
an element are satisfied ’in the mean’ through the Reissner 
variational principle. Since the state variables contain all 
quantities of interest to the designer, no further effort 
or approximation is needed to evaluate these. 

Error analysis based on the classical order of accuracy 
approach is developed. Convergence to the true solution has 
been proved and is shown by numerical experiments. The order 
of accuracy developed here is useful to obtain a ’deferred 
approach to the limit’. A few numerical results have been 
presented to show the effectiveness of the proposed model. 
Solutions have also been obtained for various shapes and 
boundary conditions of the plate and the results agree 
closely with those published in the literature. 
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CHIHJBR 7 


SHELL MIL YS IS 


7 « 1 Int rc duct ion : 

In this Chapter, finite element equations for shell 
structures will he derived through Beissner’s energy approach. 
This will be a continuation of the previous chapter. 

Many investigators have ^ready performed the stress 
analysis of shells using finite element technique with vary- 
ing degree of success. Since, many references are available 

on this study, only a fev/ recent Investigations will be nen- 

12 3 

tioned here. Clough and Johnson , Carr , Irgyris et al. and 
others have used flat , quadrilateral and triangular ele- 
ments for the shell discretization. But for curved struc- 
tures, curved elements are desirable because they eliminate 
01 * at least minimize the error involved in modelling the struc- 
ture. 

jln almost completely compatible curved triangular ele- 
ment has been demonstrated by Bonnes, Dhatt, Giroux and 
Hobicband^, Bor cylindrical shell, a completely compatible 
Curved rectangular model has been developed by Bogner, Box 
and Schmit”. Other notable contributors are Oden , Cantin 
and CloughJ. lonlinear analysis of shell structures have been 
done by Stricklin, Hsu and Pian®, Schmit, Bogner and Box^, 
Wempner^^ and Saeed Yaghmai^”*. 
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A mixed finite element model for linear as well as 
for nonlinear**^ theory of shell have been developed by Prato 
via Reissner’s principle. He has assxxmed linear displace- 
ments and moments over a triangular element. In this Chap- 
ter, the shell equations in general curvilinear coordinates 
will be discretized talcing linear distribution of all the 
fifteen unknowns as obtained through Reissner’s shell theory. 

7*2 Goveining Equation; 

Consider a curvilinear coordinate system where (x, y) 
represent axes in the plane of the shell and z normal to the 
plane, ibllowing Naghdi^^, the governing variational equa- 
tion for the shell is given by, 
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In Eq, (7*2-1), E is Yoimg’s moduliis, v is Poission’s 
ratio and h is the thickness of shell; ox and a y are Lame’s 
parameters for curvilinear surface and and Ry are princi- 
pal radii of curvatures of the middle surface; U^, Uy, Py 
and W are displacement components, whereas N , N , N , H , 

Ji j Zj Jit 

^x’ "^y* ^^xy ^yx stress resultants; p^ and Py 

are membrane loads; q and q“ are normal loads on the positive 
and negative sides of the shell element with respect to fte 
thickness of the element. Over and above, the following no- 
tations have been used in Eq. (7,2-1); 


GPi = ^ ; 0^2= 1-h2 CP^/12 Ry 

CF^ = 1 + h^ CF^/12 R^ CF^ = 1-3h^ GF^/20 Ry 

GF^ = 1 + 3h^ GF^/20R^ CFg = 8/15 - 2h^ CF^/105 Ry 

CFy = 8/15 + 2h^ CF^/105 \ 

and,H'*’ = 1+h(l/R^ + l/Ry)/2 + Ry 

ir = 1-h(l/R^ + l/Ry)/2 + h2/4R^ Ry 


(7.2-2) 
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For more elalsorate derivation of the Eq, (7.2-1), see ref- 
erence 14. 


7.3 Finite Element Derivation: 

Assume lineac distribution (in the curvilinear coordi- 
nates X, y) over a curved triangular element for all the fif- 
teen field variables. Considering a typical quantity H and 
its nodal values as , Hg and the interpolating function 

for H will be given by, 

* * * 

H = (7.3-1) 


where , 
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and are curvilinear lengths, similar to those defined 

in Eq. (6.3-4). Substituting the distribution (7.3-1) in 
Eq. (7.2-1), the following discretized form of the equations 
are obtain edj 
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7.4 Hote for Numerical Example: 

A computer program has heen developed for the shell 
aiialysis which has fifteen degrees of freedom per node. Due 
to limitation in memory space, the program can accommodate 
only a half hand width of four node points, For this severe 
restriction, it has not been possible to solve any problem 
of practical interest. However, for few special geometries, 
the resulting matrices have been checked from statical con- 
sideration. 
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ff (Ij. Ij '=®i/Sc “y^ 
ff CTj, ?j Oy) 8< 
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= If (Ti CF^ h^/12 Cy) da^/aj dx dy 
I3I = If (Ti CF^ \?/n 0 ^ ay) day/ax dx dy 
lil = If \ CFg/Eh dx dy 

i -i * * 

I34 = II 2 CFg (l 4 -v)/lh dx dy 

= // T. CFg/o^ Oy) dOy/dx dx dy 

I3I = fl F. GF^/Eh dx dy 

* * "M" 

13^ = If 2 T. GF^ (1+v)/Bh dx dy 

I3I = II (2^ “y^ ^ 

I3I = ^3 ^ 

= If 24 CF^/(l+^> )/Sb 5 dx dy 

* * -X* 

= If (Tj_ ^^4/“x “V^ '^y 

* * ”X* r? 

I4I = // 12 T^ CF^/Eh^ dx dy 

’X’ ’X" 

^3 = II 24 CF5(1+V)/E]h5 ^ ^y 

I^i = II (?4 ^^5/«x “y^ 

= // 4.5 (1 + v) T. CFg/Eb dx dy 
^4! = // 4.5 (1+v) ?. T. GF^/Eb dx dy 
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COICLUSIOH ME COMElWS ON FJTUKI! RSSSARCK 

3.1 Conclusion: 

In this sequel, following the concept forwai^ed hy 
Green and Rivlin (Ref. 4 of Chapter 1), the general govern- 
ing equations for a continuum have been derived on the basis 
of thermodynamic j.aws and material frame indifference prin- 
ciple. These equations can be directly applied for finite 
element discretization. By suitable manipulation, displace- 
ment, force 'or mixed model can be obtained separately. The 
equations essentially represent Galerkin-type approximation 
applied on Cauchy’s laws of motion and on equation of thermal 
state. Application of finite element method discretizing the 
space variables, yield, in general, a set of nonlinear ordi- 
nary differential equations in terms of the unknown field 
variables. When parametric differentiation technique is app- 
lied, the resulting modified equations become linear ordinary 
differential equations of variable coefficients in terms of 
the differentials of the field variables with respect to a 
parameter and can be conveniently integrated by any standard 
numerical technique. The unknown field variables can now be 
obxainod 'by employing any quadrature formula. In essence, 
the nonlinearity of the problem is restricted only in the 
quadrature manipulation which virtually poses no difficulty. 


The scope of application of these equations is ver^; ’.vide, 
jlctnally, v;ith minor modifications, they can he applied to any 
particular cases of continuum problems, for example, mtiLti- 
polar cases, viscoelasticity, viscoplasticity, coupled ther- 
moelasticity with dissipative properties etc. It appears that 
the applications to problems of nonlinear stability or compli- 
cated fluid flow -with thermal effects may also be possible. 
Although, in this study stresses have been expressed as the 
functionals of strains and their time rates, the constitutive 
equations expressed in the other way, i.e., expressing the 
strains aJid their rates as the functionals of stresses and 
their time rates are also feasible. This will be rather 
straight forward for mixed model, but for other cases it may 
be quite complicated. The geometrical transformations, for 
problems like plates and shells can also be done in the usual 
manner. 

In this- study, various particular cases of simpler 
types have been solved utilising this general technique. The 
results, in general, show excellent agreement for one step 
solutions, i.e., for linear ca^es. The results are also 
q-aite satisfactory for short range processes. But for long 
range solutions, the incremental step length becomes a cri- 
tical governing factor. Unfortunately, due to the restric- 
tion in available computer time, no further refinement in 
this direction has been possible. But it appears that it will 
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-v v,MC conplf'Xities than the numerical solution of 
nc'Ui iir.. lif i"' ren+ia.! equations of initial valu^^ protlcr.::. 

I.-n it i3 felt necessar:/ .xo comment on the va- 
riational. techniques commonly apt'/liel for the derivation of 
fin it-., olomt, r.t met hods. As such in the present set up, the 
variational techniques seen to he a more favourite, possible 
he cause, the discontinuities along the interelenent heun- 
daries can he easily accounted and convergency can he proved 
on the harir. of some pctrntial. Eu-^, sine: the vaimational 
technique is fundamentally an heuristic approach, the nodi- 
ficati-ons which can he incorporated there, may also be app- 
lied o:o xl;c 06 '^. of equations derived in this variant. Though, 
thiough the concept of potential it is convenient to prove 
convcrgcncu in the mean, this is not always necessary. More- 
over, for the proof, a requisite predefined norm can also he 
employed. For example, in Chapter 6 the convergency of the 
pla.tc bending mod':! has been shovn by employing the classical 
•order of accui’acy' technique, aJLtbou>~h its scope of appli- 
cability is rather limited. However, it necessitates more 
thoroub'h and extensive investigations. 

1.2 Comments on Future Research: 

There is -so much work yet to be done on continuum 
mechanics problems, that it is difficult to pent ictO. arise 
my orr- of then,. Except a few exceptions, the multipolar case 
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have remained praotioally untouohed. Moreover, the tame 
deoendent and ..nlinear e«eota on platen and nhelXB are 
also getting more and more prominanoe fer more aoc»arate 
^yeio. Aotnolly. all the prohleme dealt here aleo 

require more ir.tensiTe azid tbo rough studies. 



APPENDIX A 


A. 1 i Note on fhe Computer Programs; 

Separate computer -programs have been developed for the 
tvpcG of prcblemf? dealt in Chapters 3 to 7. Their underlying 
j-'m.c.i legt'u-e vill described here. 

Each prograro can be subdivided into the sections; 

(1) Control, (2) otiffneos, (3) Assembly, (O Elimination 
(5) Irgut/Oiitrut and (6) .Anociliaries. The execution will 
start through control. It first calls the I/O section to 
supply initial informations. Then it sets the necessary 
indices for the nodes as v.-ell ao element::, for example it ge- 
nerates the number of elementn oixrrounding a node. At the 
same time, the control calls the stiffness which supply the 
matrices through a logical tape unit. When this is over, the 
control calls the Assembly where the matrices which were 
supplied elenertfwise will be sorted nodewiae. Then the asso- 
ciated matrices are assembled for each node and immediatel;/, 
the Elirnin.allon section is called to eliir.inate the matrices 
•ppto tlie diagonru teimj malcing it an identity matrix. The left 
over portion, i.e., half band width behind the diagonal, will 
b~ -orc.ser-'/cd on a logical tape unit for the back sweep to eva- 
1.1. -r -^he iicoe.l •.inknovTis. The simple Gauss elimination prooe- 
wi;-? hot beer, rmuloyed here. During this execution lo'icpl tare 
ur,'" of T31! 7044 computer are utilised to save cieinrr:,' 
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ftice the rjodal unknowns are lonown, various other 
S"’7.i2iE7'’ factors can be calculated through subroutines which 
h'-'-.’c 'p’ 'Uped as Arocii iaries. 

This procedure may be associated with an integration 
tochr. icue courlei v.ltb s quadrature for solving the increiaeu- 
tal prciCc.'U-C. In this study, modified Euler procedure with 
trspesctdal cmlc 3 been incorporated. 

The rain disadvantage of this program is the excessive 
use O’ tapes have made it less efficient on the basis 

of executioi! time. 
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